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Behaviour of small regions of different gases carried 
in accelerated gas flows 


By GEORGE RUDINGER AnD LOWELL M. SOMERS 


Cornell Aeronautical Laboratory, Buffalo, New York 
(Received 2 April 1959) 


Small regions in a flow where the density is different from that of the surrounding 
gas do not exactly follow accelerated motions of the latter, but move faster or 
slower depending on whether their density is smaller or larger than that of the 
main flow. This behaviour cannot be quantitatively explained by treating a gas 
‘bubble’ as a hypothetical solid particle of the same density, because a gas bubble 
cannot move relative to the surrounding gas without being transformed into 
a vortex which absorbs part of the energy of the relative motion. 

To illustrate the acceleration effect, the flow velocity behind known pressure 
waves in a shock tube is compared with the observed velocity of a bubble pro- 
duced by a spark discharge. The displacement of such a bubble by a wave exceeds 
that of a flow element by more than 20 %, but the bubble density is not known. If 
the spark discharge is replaced by a small jet of another gas, a pressure wave cuts 
off a section of this jet which then represents a bubble of known density. 

A theory is developed which permits computing the response of such bubbles to 
accelerations. The ratio of the bubble velocity to the velocity of the surrounding 
gas depends on the density ratio for the two gases and on the shape of the bubble, 
but not on the acceleration. Experimental results with H,, He, and SF, bubbles 
in air, accelerated by shock waves of various strength, are presented and agree well 
with the theoretical predictions. The results apply regardless of whether accelera- 
tions are produced by pressure waves in a non-steady flow or by curvature of 
streamlines in a steady flow. Various aspects of the experimental observations 


are discussed. 


1. Introduction 

The technique of using tracers for the visualization of streamlines in a gas flow, 
or for the measurement of flow velocity, is well established. It is also known that 
such tracers do not move exactly with the flow if their density is different from 
that of the surrounding gas and if they have to undergo accelerated motions. The 
behaviour of solid particles under such conditions was described in considerable 
detail, for instance, by Wright (1951), Cady (1954), and Hoenig (1957). While such 
particles always lag in accelerated flows, gaseous particles (‘bubbles’) may lead or 
lag depending on whether their density is smaller or larger than that of the 
surrounding gas. The behaviour of gas bubbles in accelerated flows is, however, 
different from that of a hypothetical solid particle of the same density, because 


a gas bubble cannot move relative to the surrounding gas without being trans- 
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formed into a vortex which absorbs part of the energy of the relative motion. This 
effect appears to have received little attention, although its understanding is 
important not only when gas bubbles are used as tracers but also for other 
phenomena. For example, if pressure waves interact with a flame, the different 
responses of the regions of different density lead to a deformation of the flame 
front; this modifies the rate of heat release, and causes secondary waves to be 
emitted by the flame. Such waves were observed, for instance, by Markstein 
(1957), and the theoretical and experimental aspects of such interactions were 
reviewed by Rudinger (1958). 

Gaseous tracer bubbles can be produced by various methods, and such tech- 
niques, developed mainly for the measurement of constant flow velocities, were 
reviewed by Cady (1954). A method to measure not only velocities but also 
velocity distributions is due to Townend (1937), who discharged straight sparks 
between two electrodes and determined the subsequent motion and deformation 
of the spark columns. A spark column, consisting of hot gas, can be observed by 
means of a schlieren system. In an alternative method, a number of sparks, 
discharged in rapid succession, can be photographed directly. In the latter 
method, each spark tends to follow the ionized path of the preceding one, and, 
since the spark columns are carried with the flow, their displacements and 
deformations indicate the local flow velocity. This multiple-spark technique has 
also been studied by Saheki (1947), and by Bomelburg (1956), Herzog (1957), 
Herzog & Weske (1957), and Weske (1958). The errors made by neglecting the 
response of the spark column in accelerated flows are not negligible, and Herzog 
& Weske (1957) propose a semi-empirical method of correction. 

In §§ 2 and 3, shock-tube experiments are described in which the flow velocities, 
produced by pressure waves of known properties, are compared with the velocities 
of spark-heated bubbles to illustrate the large errors that may result if such 
bubbles are used as tracers in a flow. Since the effective density of a spark-column 
is not known, a new technique for producing bubbles of known density is described 
in §4, and a theory of wave interaction with such bubbles is derived in § 5. The 
theoretical and experimental results are compared and discussed in § 6. 


2. Apparatus for spark experiments 

The response of gas bubbles to accelerations is conveniently studied by 
means of a shock tube. In the first experiments, bubbles are produced by an 
electric spark discharge, and their motion is observed with the aid of a schlieren 
system. 

A vertical shock tube is used for these experiments. It has a square cross-section 
of 3in. width and a total length of about 9 ft. The diaphragm is located 3 ft. from 
the top, and the last 9in. from the bottom represent the test section which is 
provided with two glass walls for the schlieren observations. To produce shock or 
expansion waves of desired strength, the initial pressure in the test section is 
always atmospheric and the pressure in the driver section is raised or lowered as 
needed. One or more layers of photographic film are used as diaphragms and are 
punctured by a solenoid-operated plunger needle. Details of the construction of 
the shock tube have been published elsewhere (Markstein 1957). 
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A spark gap is formed between one wall of the test section and the point of 
a steel needle, introduced through the opposite wall. The spark is produced by 
charging a 1#F condenser to 400 V and then discharging it through an automobile 
spark coil in series with a 3 D 22 thyratron. The latter acts as a switch for proper 
timing control. In agreement with observations by previous investigators 
(Townend 1937, Saheki 1947 and Bomelburg 1956), reproducible straight sparks 
can be obtained if the needle electrode is negative and if the gap does not exceed 
about 2 in. To assure that the spark is discharged at right angles to the shock-tube 
wall, the latter is covered with an insulating sheet (photographic film) with a pin 
hole opposite the needle. The spark-heated column of air is clearly visible on 
schlieren photographs obtained with a short-duration spark light source. 
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FIGURE 1. Schematic arrangement of shock tube and instrumentation 
for spark experiments. 


The experimental set-up is shown schematically in figure 1. Proper timing of 
the spark discharge and of the schlieren spark is achieved by means of two 
variable-time delays. A pressure transducer is mounted near the top of the test 
section, so that the amplified signal from a pressure wave propagating down the 
tube triggers a thyratron several hundred microseconds before the wave reaches 
the spark gap. The thyratron circuit produces a single pulse which, in turn, triggers 
the spark discharge after a desired delay. A pick-up coil, consisting of a few turns 
of wire wrapped around the high-tension lead of the spark, produces a signal 
which, through a second thyratron pulse generator and the second delay unit, 
triggers the schlieren spark. 


3. Measurements of flow velocity behind pressure waves 

In the first series of experiments, the displacement of the spark column in the 
wake of a weak shock wave is measured as a function of time and compared with 
the corresponding displacement of an element in the surrounding gas. The spark 
is discharged shortly before the arrival of the wave at the spark gap, and aschlieren 
photograph is taken after some interval. Figure 2a, plate 1, shows a typical 
record. 

Such photographs can be evaluated with reference to the wave diagram in 
figure 3a which shows the motion of the shock wave, the spark column, and 

11-2 
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a corresponding gas element. Let x, be the distance of the shock wave from the 
spark gap at the time the schlieren photograph is taken. If this distance is plotted 
as a function of the delay setting for the schlieren spark, one obtains a straight line 
which intercepts the line xz, = 0 at the instant that the shock wave passes the 
spark gap. Let At be the time interval in which the wave has travelled the distance 
xy. These data determine the shock velocity and, together with the initial speed of 
sound in the test section a,, the shock Mach number W,. For the weak shocks used 
in these experiments, the results agree well with calculated data based on the 
theory of an ideal shock tube. 
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FIGURE 3. Wave diagrams for spark experiments with the spark discharged: 
(a) before, and (b) after, passage of a shock wave. 


The spark bubble remains at rest until the shock wave arrives; it then travels 
the distance 2, during the time At. During the same time, a gas element that 
initially was located at the spark gap travels the distance z,. If the flow were 
accurately represented by the bubble motion, the ratio of the two displacements 
should be unity. Since At = xg/ug = xg/(4y Mg), where uw, is the flow velocity of 
the gas, the actual displacement ratio is given by 


(1) 
The first factor on the right-hand side of this equation can be obtained directly 


from the schlieren photographs, and the second factor varies only with Mg, 
since wu, is related to M, and a, by 


. (a1 : (2 
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where y is the ratio of the specific heats. The displacement ratios observed for 
shock Mach numbers of about 1-12, 1-22 and 1-28, corresponding to shock-pressure 
ratios of 1-29, 1-58 and 1-74, respectively, are plotted in figure 4. This plot indi- 
cates that the velocity of a spark-heated bubble immediately after interaction 
with a shock wave is about 60 % higher than the flow velocity. Subsequently, the 
bubble velocity decreases and gradually approaches the flow velocity. In spite of 
the scatter of the data, the uncertainty of the measured bubble displacements is 
considerably smaller than the deviation of the displacement ratio from unity, as 
can be seen in figure 2a, plate 1, where the value of x, is also indicated. 
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Time, At, msec 
FicurE 4. Displacement ratios of spark bubbles for shock waves of various strength as 
a function of the travel time of the bubbles. Spark is discharged either before (open 
symbols) or after (filled symbols) passage of the shock wave. Shock pressure ratio: ©, 1-29; 
A; @, 1-58; ©), 1-74. 


To verify that the observed deviations of the displacement ratio from unity are 
caused by the interaction of the shock wave with the spark volume, the bubble 
displacement is also determined in experiments in which the spark is discharged 
after the shock wave has passed the spark gap. A typical schlieren photograph is 
shown in figure 25, plate 1, and the wave diagram corresponding to these experi- 
ments is given in figure 3b. The travel time of the bubble At is determined by the 
setting of the second delay unit (see figure 1), but a correction must be applied 
because of unavoidable fixed delays in the spark circuits. This correction is 
obtained by plotting x, as a function of the delay setting and extrapolating to 
Xp = 0. The displacement ratio x,/x, is then obtained from the measured value 
of xp and from xq = u,,At. The results are indicated in figure 4 by the filled-in 
symbols and verify that the bubble moves with the flow as long as accelerations 
are absent. 

From the data presented in figure 4, it appears that the excess bubble velocity 
depends only little, if at all, on the strength of the shock wave. Although the 
shock accelerations are extremely high, the shock Mach number is varied over 
a small range only. The effect of accelerations of a different order of magnitude 
can be studied by experiments in which the spark column is accelerated by an 
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expansion wave. Such waves are produced by partially evacuating the driver 
section of the shock tube, and their properties can be determined from oscilloscope 
records of the pressure variations at the location of the pressure transducer (see 
figure 1). The expansion wave is reflected from the closed end of the shock tube 
and leaves a region of reduced density and zero velocity. If the spark is discharged 
before arrival of the wave and the schlieren photograph taken after passage of the 
reflected wave, the bubble displacement x, is readily compared with the corre- 
sponding displacement of a gas element x,. To test for proper timing of the events, 
the starting pulse from the second time-delay (see figure 1) is also used to produce 
a signal on a second channel of the oscilloscope to mark the instant of the spark 
discharge. The instant of the schlieren spark is similarly indicated on the pressure 
record. Figure 5, plate 1, shows a schlieren photograph with the corresponding 
oscilloscope record, where it can be seen that the spark is discharged before the 
pressure starts to decrease, and that the photograph is taken after the pressure 
again has become steady. The interval between the spark discharge and the 
taking of the photograph is about 4-5msec. A spark column of this ‘age’ no 
longer appears as a straight line but as a bubble of irregular shape. The mean of 
the maximum and minimum distances of the bubble boundaries from the spark 
gap is, therefore, taken as xp. 

During the passage of the incident and reflected expansion waves, the density 
decreases from py to p,,, and this change can be computed from the corresponding 
pressure change since the flow is isentropic. Let L be the distance of the spark gap 
from the closed end of the shock tube. The value of x, is then related to the 
densities and pressures by 


tgth _ Po _ (Po ) . 
L Px Dol 


(3) 


The pressure ratio p,,/p) for these experiments is in the neighbourhood of 0-7. 
For the record in figure 5, plate 1, p,,/p) = 0°66, and the value of x, computed 
from equation (3) is indicated in the figure. Because of the uncertainty of x, the 
displacement ratios scatter considerably, and the observed values range from 
about 1-15 to 1-25. Since gradual cooling of the bubble during its acceleration and 
deceleration by the waves tends to reduce the excess displacement, these values 
seem sufficiently close to those observed for shock accelerations (figure 4) to 
support the previous conjecture that the displacement ratio is independent of the 
applied acceleration. 


4, Experiments with bubbles of known density 

The experiments described above give a good indication of how pressure waves 
interact with bubbles the density of which is different from that of the surrounding 
gas, but the experimental results cannot be used to check a theoretical analysis 
because the density of a spark-heated bubble is not known. A technique was 
developed, therefore, to produce bubbles of known density. 

At first, soap bubbles filled with various gases were tried, but this approach was 
soon abandoned. The displacement of large bubbles amounts to only a few bubble 


diameters and cannot be measured accurately because of changes of bubble shape. 
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For small bubbles, the mass of the soap film represents a significant fraction of the 
total mass, so that the bubble density becomes uncertain. 

In the technique finally adopted, a fine jet of some gas, discharged across the 
shock tube, represents a ‘bubble’ of known density. For this purpose, the needle 
electrode of the spark gap is replaced by a piece of hypodermic tubing of 0-033 in. 
internal diameter. A second tube of about 0-25in. diameter is mounted in the 
opposite wall of the test section to receive the jet discharged by the fine tube. The 
larger tube is connected to an aspirator pump to remove the jet flow before it can 
contaminate the air in the test section. The gap between the two tubes is about 
0-3in. so that a passing shock wave cuts off this exposed length of the jet. The 
motion of this jet element can then be observed with the aid of a schlieren system. 

Figure 6, plate 2, shows four schlieren photographs obtained by this technique 
for a shock-pressure ratio of about 1-3. The gases used are hydrogen, helium, and 
sulphurhexafluoride. Such records can be evaluated by the procedure already 
described in connexion with the spark experiments, and the values of 2,, corre- 
sponding to the experimental conditions are marked in the figure. It is seen that 
the displacements of the bubbles are larger than 2, for the light gases H, and He, 
and smaller for the heavy gas SF,. 

One interesting feature of these photographs is that all bubbles are inclined 
although the jet is discharged horizontally. This phenomenon is readily explained. 
The pure gas that leaves the jet tube mixes with air as it advances across the test 
section. Consequently, the effective density of the bubble changes continuously 
from that of the pure gas to that of a mixture of this gas with air. The bubble end 
nearest the jet tube, therefore, travels faster than the other end for H, and He, and 
slower for SF,. Because of the rapid deterioration of the SF, bubble, which is 
discussed in § 6, the inclination of this bubble is clearly visible only on photographs 
taken early enough after passage of the shock wave (lower right picture of 
figure 6, plate 2). 

The extent to which the jet mixes with air, and therefore the inclination of the 
bubble, depends on the diffusion rate and on the jet velocity which controls the 
time available for mixing. No attempt was made to keep the jet velocity the same 
for all gases, and the photographs indicate that the velocity of the H, jet was 
considerably higher than that of the He jet, because the He bubble is more 
inclined than the H, bubble in spite of the 10% larger diffusion coefficient for 
H,-air as compared to He-air (Hirschfelder, Curtis & Bird 1954). 

In view of the foregoing discussion, all bubble displacements must be measured 
on the side of the bubble nearest to the jet tube. Experimental results obtained 
with the described technique are discussed in § 6. 


5. Theory of wave interaction with gas bubbles 

For an observer moving with the gas, the acceleration by a wave is equivalent 
to that of a gravitational field. In this field, two forces act on a bubble in the gas, 
namely, its weight in the equivalent gravitational field and a buoyancy force equal 
to the corresponding weight of the gas displaced by the bubble. Those two forces 
act in opposite directions and are of equal magnitude only if the density of the 
bubble is the same as that of the surrounding gas. If the bubble moves because the 
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resultant force is not zero, its surface becomes a sheet of vorticity. Subsequently, 
this sheet rolls up with the result that the entire bubble is transformed into 
a vortex. The effect of this transformation corresponds to a drag on the transla- 
tional bubble motion relative to the gas which must supply the energy absorbed 
by the vortex flow. 

Vortex formation is not readily apparent in the described experiments because 
the bubbles are too small to reveal sufficient details of the flow, although figures 
2a and b, plate 1, indicate the more ‘turbulent’ character of a bubble that has 
undergone interaction with a shock wave as compared to the smooth surface of a 
bubble that has not interacted. The development of a vortex can be clearly seen in 
figure 7, plate 3, which shows three stages of the interaction of a shock wave with 
a flame bubble. These schlieren photographs are taken from a series obtained by 
Markstein (1958) in the course of another investigation. A combustible mixture 
in the test section of the shock tube is ignited by a spark between the electrodes, 
visible at the centre of each photograph, and a bubble of burned gas is contained 
within the flame front. The latter advances so slowly that its progress during the 
time of interest here is unimportant. The first photograph shows the approximate 
spherical flame bubble just before interaction with a downward-travelling shock 
wave. The second photograph shows the wave emerging from the bubble and also 
a reflected wave. In the third photograph, the waves have moved away from the 
bubble which has now, clearly, been transformed into a vortex ring. 

The following analysis by-passes a detailed study of this vortex formation and 
yields the final translational velocity of the bubble directly on the basis of a few 
plausible assumptions. Consider, first, accelerations by shock waves. The time 
during which a shock wave passes through the bubble, i.e. the time during which 
acceleration takes place, is so short that the vortex must essentially be formed in 
a second, subsequent stage. This view is supported by figure 7, plate 3, where no 
indication of vortex formation can be seen at the time when the shock wave has 
left the bubble. The analysis will, therefore, proceed in two steps, corresponding 
to the initial acceleration of the bubble and to its subsequent transformation into 
a vortex. Later, it will be shown that, within the accuracy of the assumptions 
made, the same two-step analysis may also be applied to waves, other than shock 
waves, for which acceleration and vortex formation take place simultaneously. 

During the acceleration stage, the bubble distortion may be ignored (centre of 
figure 7, plate 3), and it seems permissible to treat the initial phase of the motion 
in the following way. Let the subscript P identify the quantities associated with 
the bubble motion immediately following the passage of the shock, and let the 
subscript @ identify the quantities associated with the ambient gas particles. 
Assuming that one may neglect the modification of the pressure field caused by 
the presence of the bubble, the impulse per unit volume J which the bubble 
experiences during the shock passage is precisely that which the ambient gas 
would have experienced, namely, J = pgug, where wu,; is the gas velocity after the 
shock has passed. 

This impulse must accelerate the bubble from rest to the speed wp, and must 
also accelerate some of the surrounding gas from the speed wu, to a velocity 
distribution associated with the motion of the bubble. Since the latter is not 
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distorted during this early phase, the bubble may be treated as if it were a solid 
particle. One may thus invoke the well-known concept of ‘apparent additional 
mass’ (virtual mass) for such a situation and obtain 


I = pgtg = ppuptkpg(up— ug), (4) 


with k as the apparent additional mass fraction (i.e. inertia coefficient). For 
aspherical bubble, & = 0-5, and for a long circular cylinder moving at right angles 
to its axis, k = 1-0 (Lamb 1945, p. 155). Equation (4) yields directly 


u obits u (5) 

= ty vo 

 @+k © 

where o = PP PB (6) 
Pa Pa 


is the density ratio. For real solid particles, such as are customarily used for flow 
visualization, o is of the order of 1000, so that w, is practically zero. For @ less 
than unity, the particle velocity exceeds that of the gas. The foregoing part of the 
analysis is closely related to calculations described by Birkhoff & Caywood (1949) 
who considered accelerated motions of gas bubbles in a liquid. 

The second step of the analysis deals with the transformation of a gas bubble 
into a vortex as a result of the initial motion relative to the surrounding gas. Since 
the kinetic energy of the initial motion must supply the energy absorbed by the 
vortex, the relative velocity is reduced from up—wug to Up—Ug, where up is the 
desired final bubble velocity. Let the ratio of these relative velocities be denoted 


by B, so that Up—Ug = B(up- Ug). (7) 


The quantity £2 thus represents the fraction of the kinetic energy of the initial 
relative bubble motion that is available for the final relative translational motion. 
Substitution of equation (5) into equation (7) yields the final bubble velocity 
relative to the duct as il 
bn = (+2 Ue. (8) 
where / is yet to be determined. 

Figure 7, plate 3, indicates that a spherical bubble is transformed into a vortex 
ring, and it seems plausible that an infinitely long cylindrical bubble would 
become a linear vortex pair. These two cases can be treated theoretically, and the 
bubble motions observed in the described experiments should then fall between 
those corresponding to these extremes. 

The calculations for a vortex ring are based on a study by Taylor (1953). He 
considers a flat circular disk that is suddenly set in motion along the axis of 
symmetry, in a perfect fluid, and without separation of the flow. Ifthe material of 
the disk were then suddenly ‘dissolved away’, the flow could also be considered 
as being caused by a distribution of vortex lines over the plane of the disk. The 
strongest vortex lines near the edge of the disk tend to roll up the weaker ones 
until, finally, a vortex ring remains for which it is assumed that the vorticity is 
uniformly distributed over a core of small diameter compared to the diameter of 
the ring. Such a vortex cannot remain at rest in the surrounding gas, but moves 
with a velocity that depends on its dimensions and strength. To apply Taylor’s 
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results to the present problem, the assumption is made that this vortex ring 
represents a reasonable approximation of the vortex into which a spherical gas 
bubble is transformed, provided the initial velocity of the bubble is equal to that 
of the disk. Taylor derives the ratio of the velocity of the vortex relative to the 
fluid to the initial velocity of the disk. This ratio has been defined here as /, and 
Taylor’s results is # = 0-436. 

If the circular disk in Taylor’s analysis is replaced by an infinite lamina, 
a linear vortex pair is formed which, by analogy with the preceding case, is con- 
sidered to be an approximation to the vortex into which a cylindrical bubble is 
transformed. This analysis is presented in the Appendix and yields # = 0-203. 

The foregoing treatment gives no indication of the time required for the vortex 
to form, but figure 7, plate 3, shows that the vortex is fully developed before the 
bubble has travelled one bubble diameter. Since the displacements of interest 
have a magnitude of many such diameters (see figure 6, plate 2), it may be 
assumed that the vortex motion is instantaneously established. As a consequence 
of this assumption, equation (8) may be applied also to waves other than shock 
waves. Such waves can always be approximated as closely as desired by a sequence 
of discontinuous wave elements for each of which the ratio of the velocity 
increments Aw,,/Au,, has the same value. If the two steps of the analysis were 
applied to each element consecutively, the final result would thus be the same as if 
the analysis had been applied to the entire wave, provided only that the bubble 
density is not significantly changed by diffusion during the time of wave inter- 
action with the bubble. 

Two possible sources of error in the analysis should be discussed. Buoyancy of 
a bubble in the gravitational field of the earth causes a bubble motion that is 
superposed on the motion produced by a wave. Although the acceleration of this 
motion is several orders of magnitude smaller than the wave acceleration, it acts 
during the entire experiment while the wave acts only during its passage across 
the bubble. For the conditions of the present experiments, it is readily verified 
that this effect is utterly negligible if the travel time of the bubble is less than 
Il msec, and even for travel times of 5msec, as in the experiments with the 
expansion waves (figure 5, plate 1), the displacement caused by gravity alone is 
less than 3° of the total bubble displacement. 

The other effect to be considered is related to differences in compressibility of 
the two gases involved. In the derivation of equation (8), it has been implied that 
the density ratio remains constant during the acceleration of the bubble. This 
condition is satisfied only if the gases have the same compressibility. If the 
compressibilities are not equal, the density ratio changes during passage of the 
wave from o, to a, which is given by 


(9) 


where p,,/po is the pressure ratio across the wave. The relation for isentropic 
changes of state may be used because only weak waves are considered here. Since 
an effective value of the density ratio must lie between o, and c,, the influence of 
this uncertainty of ¢ on uw, in equation (8) is less than 0-5 % for the conditions of 
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the reported experiments (p/p, < 2, Yg = 1:40 for air, yg = 1-10 for SF,, and 
Vp = 1-67 for He). 

In the described experiments, the ratio of bubble length to its diameter is 
about 9. If this shape is considered as an elongated ellipsoid of revolution, the 
inertia coefficient (Lamb 1945, p. 155) is approximately 0-94. In view of the 
uncertainty of the effective shape of the bubbles, the value k = 1-0, corre- 
sponding to an infinite cylinder, is used. The effect on w, of such a variation of k 
in equation (8) is unimportant. 

















UB/UG 
poo a 
Gas o Up/Ug Vortex ring Vortex pair 
H, 0-0695 1-87 1-38 1-18 
He 0-138 1-76 1-33 1-15 
SF, 5-03 0°33 0-71 0-86 
TABLE 1. Theoretical results 
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Figure 8. Observed and theoretical displacement ratios for the interaction of shock waves 


with gas bubbles of different densities (see also Table 1). —, Vortex ring; . vortex 
pair. Shock pressure ratio: ©, 1:29; A, 1-58; [), 1-74. 


Experimental results should lie between the values predicted by equation (8) 
for the two values of /, corresponding to a vortex ring and a linear vortex pair, 
respectively. For the gases used, the initial and final velocity ratios, wp/ug and 
p/Ug, as computed from equations (5) and (8), are collected in table 1. 

Figure 8 shows the theoretical and experimental displacement ratios plotted 
as a function of the time elapsed after interaction of the shock wave with the 
bubble. The pressure ratios of the shock waves used in the experiments are 1-29, 
1:58 and 1-74, and the experimental points are identified accordingly. The 
theoretical values are independent of time and are thus represented by horizontal 


lines. 
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6. Discussion 


The effect of accelerations in a gas flow on the motion of gas bubbles of different 
density has been studied both theoretically and experimentally. In the present 
investigation only one-dimensional flows accelerated by pressure wave are con- 
sidered, but it should be emphasized that it is irrelevant whether the accelerations 
occur in a non-steady flow or are caused by curvature of streamlines in a steady 
flow. The behaviour of gas bubbles in an accelerated flow is quite different from 
that of solid particles and cannot be explained by the difference in density alone. 
If a shock wave passes over a solid particle, which may be about 1000 times as 
dense as the surrounding gas, the initial particle velocity is practically zero. 
Subsequently, the particle speeds up under the influence of its aerodynamic drag 
until it moves with the flow. A gas bubble, on the other hand, acquires almost 
instantaneously, a velocity that may be larger or smaller than the gas velocity, 
depending on whether the density ratio is smaller or larger than unity. This 
velocity is, however, different from that of a hypothetical solid particle of 
the same density because a gas bubble in relative motion with respect to its 
surroundings is transformed into a vortex which absorbs part of the energy. 
Subsequently, the relative bubble velocity gradually decreases as a result of 
viscous dissipation. 

The theory is based on the assumption that a spherical bubble is transformed 
into a vortex ring, and a cylindrical bubble into a linear vortex pair. One obtains 
for these limiting cases bubble velocities that differ by less than 20 °% for the gases 
used here, and the velocities for other bubble shapes should fall between these 
rather narrow limits. The experimental data in figure 8 are in good agreement 
with this requirement and thus support the main features of the theory, namely, 
that the vortex is formed practically instantaneously (before the bubble has 
travelled one bubble diameter), and that the velocity and displacement ratios are 
independent of the magnitude of the acceleration. 

According to equation (7), #? represents the fraction of the energy of the initial 
relative bubble motion that is available for translational motion after a vortex is 
formed. The effective value of # should lie between 0-203 and 0-436, so that the 
major fraction of the energy is absorbed by the vortex, and only about 10% 
remain for the translational motion. 

It is interesting to compare the results for gas bubbles with those for solid tracer 
particles, for example, those given by Hoenig (1957). A shock wave of Mach 
number 10, propagating in argon of 273°K, at a pressure of lem of mercury, 
produces a flow velocity of 2-42 x 10° cm/sec. A particle of 10-*cm diameter and 
a density of 3 g/cm? remains practically at rest immediately after passage of the 
shock wave, and reaches 95 °% of the gas velocity after about 0-15 msec. After this 
time, the particle follows the flow closely. By comparison, the decrease of the 
relative velocity of a gas bubble behind a shock wave is so slow that it is barely 
noticeable in figure 8. The displacement ratios for H, and He show a tendency to 
decrease and those of SF, to increase with time, but the decay of the relative 
motion is still smaller than the scatter of the experimental points even after 


0-3 msec. 
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The displacement of a gas bubble that initially is moving with the flow can 
always be considered as the vector sum of two displacements, one that would be 
observed if the flow velocity remained constant, and the other caused by accelera- 
tions. The former exactly represents the corresponding motion of the flow, while 
the latter differs from that of the flow by the factor 1 + #(1—0)/(a +k), according 
to equation (8). This factor is independent of the manner in which the acceleration 
takes place as long as o may be treated as a constant. If the acceleration lasts long 
enough for o to change significantly as a result of diffusion and/or heat conduction, 
the foregoing constant factor no longer applies. Dissipation of vorticity because 
of viscosity also reduces the relative bubble motion. The effect of diffusion, which 
causes the sloping appearance of the bubbles in figure 6 (plate 2), has been dis- 
cussed in §4. Although one is dealing there with variations of o in space rather 
than in time, these effects are closely related. 

Another interesting feature of the behaviour of gas bubbles is the instability of 
bubbles that are heavier than their surroundings. Figure 6, plate 2, illustrates the 
rapid disintegration of SF, bubbles, while H, and He bubbles maintain their 
appearance for considerably longer times. Similar observations about the stability 
of buoyant vortices were reported by Turner (1957), who also provided the ex- 
planation that a light vortex core stabilizes the motion in a manner comparable to 
that of a stable stratification of two gases in a gravitational field, except that, in 
the vortex, gravity is replaced by the centrifugal force due to rotation. 

Figure 4 shows that the displacement ratio of spark-heated bubbles decreases 
rather rapidly, while figure 8 indicates almost constant values for bubbles of the 
three gases used. The main difference between these experiments is that spark- 
heated bubbles are much hotter than the surrounding gas, while the temperature 
of the gas bubbles differs only a little if at all. It is readily verified that the 
compression produced by the strongest shock waves used in the experiments 
causes differences between bubble and gas temperatures of 0 °K for H,, 23 °K for 
He, and — 23 °K for SF, if the initial air temperature is 300 °K. Heat conduction 
and viscosity effects should, therefore, become significant mainly for the spark- 
heated bubbles. The consequences of these effects are not included in the theory, 
and only a tentative explanation of the observed behaviour is offered here. It 
should first be noted that, although the density of a spark-heated bubble is not 
known, the value of the displacement ratio computed from equation (8) cannot 
exceed 1+, corresponding to o = 0, and for k = 1-0. Extrapolation of the 
experimental data in figure 4 to the origin of the time scale should, therefore, 
yield a value between about 1-2 and 1-4, depending on the effective value of #, but 
the experimental data indicate an initial displacement ratio in the neighbourhood 
of 1-6. The observed results for small values of At seem, therefore, too high to 
represent the decay of a vortex. The data lie between the values corresponding to 
a fully developed vortex, as given by equation (8) and those corresponding to the 
initial bubble velocity before a vortex is formed, as given by equation (5), but 
there seems to be no reason why the transformation of a hot bubble into a vortex 
should be slower than that of a cold bubble. The high bubble velocities may be 
explained by considering the effect of cooling on the vortex. The propagation 
velocity of a vortex relative to the surrounding gas is proportional to the circula- 
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tion and inversely proportional to the size of the vortex (see, for example, 
equation (A5) of the Appendix). While cooling does not affect the circulation, it 
vauses the vortex to contract and, therefore, to speed up. According to this 
hypothesis, the behaviour of spark-heated bubbles, as shown in figure 4, is an 
indication of the extent to which the slowing down of the bubble because of vortex 
iormation is counteracted by a speeding up due to cooling. 

The large viscosity of air at elevated temperatures should cause a rapid dissipa- 
tion of vorticity. An analysis of the slope of the relation between the displacement 
ratio and time (figure 4) shows indeed that after about 0-2 msec, the bubble 
velocity is already reduced to within a few per cent. of the flow velocity. It appears 
therefore that, immediately after interaction, the behaviour of a spark-heated 
bubble is significantly affected by cooling, while later, the viscosity effects 
dominate. 

If spark-heated bubbles are used for tagging a flow, a correction for acceleration 
effects may have to be made to avoid appreciable errors, particularly immediately 
after accelerations. For bubbles formed by different gases, the acceleration effects 
can be greatly reduced by selecting density ratios as close to unity as the sensi- 
tivity of the detection system permits. In the described experiments, rather 
extreme conditions are used to demonstrate the phenomena clearly and to 
provide a sensitive test for the theory. The described results not only have 
a direct bearing on the use of tracer bubbles in gas flows, but they may also 
contribute to the understanding of various phenomena in flows of non-uniform 
density. 


The authors wish to thank their colleagues G. H. Markstein and A. Gail for 
stimulating discussions in the course of this work. This work was started as an 
internal research project of Cornell Aeronautical Laboratory, and continued 
under the auspices of Project SQUID which is supported by the Office of Naval 

Xesearch under Contract Nonr-1825(25). 


Appendix 
Impulsive motion of an infinite lamina and subsequent formation 
of a linear vortex pair 
Taylor (1953) analysed the formation of a ring vortex under the condition that an 
impulse is given to a circular disk and the latter is then ‘dissolved away’. 
A similar analysis is carried out in the following for an infinite lamina. 

Assume that a lamina of width 26 is suddenly moved with a velocity U at right 
angles to its plane and that no flow separation occurs. A discussion of the 
resulting potential flow is given by Lamb (1945, p. 85), and one obtains for the 
kinetic energy of the motion 7H _ Lpb? U2, (Al) 
and for the potential of the flow at the surface of the lamina 

= 1(b2 — 423 2 
p= +U(b—-y?)?, (A2) 


where y is the co-ordinate in the plane of the lamina, measured from its centre. The 


impulse of the motion is given by 


P =dT/dU = 7pb?U, (A3) 
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and the surface velocity by 
dg/dy = + Uy(b?—y?)-3. (A4) 


If the lamina were suddenly dissolved away, the flow could be considered as 
being produced by a distribution of rectilinear vortices over the surface of the 
lamina. These vortices tend to roll up to form a linear vortex pair with the centres 
of the cores located near the edges of the lamina where the vorticity is strongest. 
It is assumed that the vorticity is uniformly distributed over the cores. Let the 
distance between the centres of the cores be denoted by 2B. The velocity of 
a vortex pair relative to the surrounding gas, and its impulse, are given by 


amb 1945, pp. 221, 229) 2 
(Lamb 1945, pp. 221 ) a (A5) 
and P’ = 2BKp, (A6) 


respectively, where the circulation around the cores is + K and — K. With the aid 
of equation (A 4), one obtains for the circulation 
*b - 
K=2| Uy(b?—y?)-2dy = 2Ub. (A7) 
/0 
Xolling up of the vortex elements does not affect the impulse of the motion, so 
that P’ = P and equations (A3), (A6) and (A7) yield 


B = inb. (A8) 
If this result is substituted into equation (A 5), one obtains 
vy = (2/7?) U = 0-203 U. (A9) 


In view of the hypothesis made for the motion of a gas bubble, U and v must be 
identified with the velocities of the initial and final bubble motion relative to the 
surrounding gas, so that U = up—u,andv = ug—Uug. Equation (A9) thus yields 
the final result # = 0-203. 
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(b) 
Figure 2. Typical schlieren photographs of spark columns with the spark discharged 
(a) before, and (b) after, passage of a shock wave (pressure ratio 1-58). Arrows indicate 
calculated displacements of an unheated air element that initially was at the location of the 


spark gap. 
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Figure 5. Schlieren photograph of spark bubble after acceleration by expansion wave and 
its reflexion from the closed bottom of the shock tube (overall pressure ratio 0-66). Pressure 
record is taken a few inches above spark gap. Arrow indicates calculated displacement 
of an unheated air element that initially was at the location of the spark gap. 
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FIGURE 6. Typical gas jet /shock-wave interactions 
(shock pressure ratio 1-58). 
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Figure 7. Interaction of a shock wave (pressure ratio about 1-3) with a bubble 
of burned gas (stoichiometric butane air mixture). 
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Taylor instability of finite surface waves 
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The instability of the accelerated interface between a liquid (methanol or carbon 
tetrachloride) and air has been investigated experimentally for approximate 
sinusoidal disturbances of wave-number range from well below to well above the 
cut-off. The growth rates are measured and compared with theoretical results. 
A third-order theory shows the phenomena of overstability which is found in the 
experimental results. Some measurements of later stages of growth agree 
moderately well with the available theory and disclose some additional pheno- 
mena of bubble competition, Helmholtz instability with transition to turbulence, 
and jet instability with production of drops. 


1, Introduction 

If a glass were held upside down with water held magically, in its usual position, 
and then released, the water would fall out. This process occurs in underwater 
explosions, in the acceleration by combustion of liquid fuels, in the acceleration 
of stellar material, and, perhaps, in the molten layer, if any, on the leading edge 
of a blunt re-entry body. These acceleration-unstable liquid surface phenomena 
were studied by G. I. Taylor in the early 1940’s with the first publication in 1950. 
He showed that in the absence of viscosity and surface tension all small surface 
irregularities are unstable and grow in time. 

An initial test of the correctness of Taylor’s ideas, apart from comparison with 
explosion photographs where the accelerations were uncertain, was made by 
Lewis (1950). As discussed later, the agreement was well within experimental 
scatter, but the investigation was limited to a wavelength (or acceleration) range 
in which the effect of surface tension was slight. 

That surface tension and viscosity have modifying effects is clear, and it was 
shown by Bellman & Pennington (1954) that surface tension produces a cut-off 
wave-number (or wavelength) above which the surface disturbance merely 
oscillates and below which the disturbances grow. Viscosity, on the other hand, 
merely decreases the rate of growth at low wave-number and causes damping 
of the oscillatory solutions at high wave-numbers. 

Allred & Blount (1953) with apparatus similar to that used by Lewis have, by 
using two fluids with densities close to each other to reduce the growth rate, 
studied the stability in a range where surface tension would be expected to play 
a part. Their results indicated a growth rate appreciably lower than the rate 
indicated by Bellman & Pennington’s theory. 

!xperimental results show that the linear growth discussed in the theories can 
12 Fluid Mech. 7 
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apply only for very short times. The liquid surface rapidly changes to liquid 
spikes falling into the gas and gas bubbles rising* through the liquid. 

Attempts to analyse these later developments have been made by Ingraham 
(1954) analytically, and Pennington et al. (1953) numerically. Both analyses 
showed the type of distortion of sinusoidal waves as observed, but neither carried 
the work to very large displacements. 

Certain large displacement effects have been examined. Davies & Taylor 
(1950) studied the constant rate of rise of bubbles through a stationary liquid, 
while Birkhoff (1954) and Garabedian (1957) have computed with high precision 
the bubble-rise velocity. 

On the other hand, the spikes have been shown analytically by Carrier (1953) 
and numerically by Pennington et al. (1953) to be simply left behind. Thus 
relative to the acceleration field the spikes are simply in free fall. 

These very approximate analyses fall far short of describing the later stages of 
development. The régime established after spikes and rising bubbles have 
appeared could remain until the spikes ‘hit the bottom’, while the bubbles ‘rise 
to the top’; but in fact such a motion is found unstable. The rising bubbles com- 
pete with one another, the large ones growing at the expense of the small ones. 
Some indications of this type of behaviour is given in figure 1 (plate 1). The gas- 
liquid interface on the sides of the spikes and bubbles are in relative (tangential) 
motion which produces additional surface instability of the Helmholtz type. 
These experimentally observed phases have been recognized by Lewis (1950) and 
Birkhoff (1954), but no detailed attempt to obtain analytical nor experimental 
knowledge of the processes has been made. In particular the final stage of mixing 
could only be of interest in its various statistical properties (it being hopelessly 
complex for detailed description), yet it could be analysed only by a statistic 
which was strongly biased by spacial and time variations of very large magnitude. 


2. Growth of surface waves 
A. Theoretical consideration—finite waves with surface tension 


For the purposes of analysis, we shall treat the medium as a non-viscous, incom- 
pressible fluid. Surface tension is included. The upper fluid is taken as air whose 
density we may neglect relative to that of the liquid. The flow-field is assumed to 
be two-dimensional, with the 2’-axis taken in the unperturbed plane of the inter- 
face. For simplicity, we shall assume that the motion of the whole system is 
started from rest and the initial surface disturbance is a simple sinusoidal standing 
wave with amplitude 7, and wave-number k’. For the proper interpretation of the 
results, all the physical parameters concerned will be expressed in dimensionless 
forms. As standards of measurement we shall take the wave-number k’ and the 
wave velocity ,/(9'/k’), of the initial disturbances, as scales of length and velocity, 
respectively, where g’ is the acceleration relative to a force-free frame of reference. 
pq’ is a virtual gravity force taken as positive in a direction normal to the surface 


* The terms ‘rising’ and ‘falling’ refer of course to the familiar gravity effects on an 
upside down glass of water. For the general case, ‘downwards’ is the direction of the effective 
gravity component, directed for instability from the heavier toward the lighter medium. 
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from the heavier toward the lighter medium. Thus for liquids at rest on the 
surface of the earth, g’ = — 981 cm/sec?. 
The dimensionless independent variables of the problem will be taken as 


saz, y=, t=(+9F1). (2.1) 


In terms of these dimensionless variables, the governing equations and the 
boundary conditions for the velocity potential d(x, y, t) (defined by u = — grad ¢) 
and the surface displacement 7(~,t) in the y-direction are: 
Governing equations: 

Pex + Pyy =0 


for oO>r>-O, Wo>y>-w, t20. (2.2) 


Boundary conditions: 
The liquid surface moves with the liquid, 


Nt Ne Pet Py = 0. (2.3) 
The pressure in the liquid differs from atmospheric by the effect of surface 
tension and the curvature, 


—19— $+ 482+ 93) = [tae + 92)-*) (2.4) 

on y = (x,t). 
Initial conditions: n(x, 0) = Ny cose, (2.5) 
n(x, 0) = 0. (2.6) 


In equation (2.4) k is given by, 


where k;, the cut-off wave-number, is given by 


, 0g’ — 
ke = If ; (2.8) 


where 7' is the surface tension. The parameter k can also be interpreted as the 
ratio of the surface tension forces to the gravity forces. 

In order to proceed analytically, let us assume that the surface disturbance 
admits an expansion of the following type: 


ie @) 


n= < N0N(x, t). (2.9) 


- 
If y, and all its derivatives are further assumed to be of the same order of 
magnitude, it follows from equation (2.3) that 
i 8) 
= L ,- (2.10) 
r=1 
Since equation (2.2) is linear, it must be satisfied by every ¢,. Our general 
scheme of solution is to introduce equations (2.9) and (2.10) into (2.3), (2.4), (2.5) 
and (2.6) and equate like powers of 7, to find the corresponding boundary and 


initial conditions for ¢, and 7,. The values of ¢, and its derivatives required at 
12-2 
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y = (x,t) are found from the values at y = 0 by use of Taylor expansions. We 
then determine each ¢, and 7, in succession to complete the solution. First, 


we have 
[oosx for r= 1,) 


n(x, 0 2.11 
ne, 0) Lo for r+1.J “— 
The following notation will be used throughout: 
Net = st for all r. (2.12) 
a 


Then, the boundary conditions obtained by the above process from equa- 
tions (2.3) and (2.4), wherein all quantities are evaluated at y = 0, are: 


MattPr,y = 0, (2.134) 

I+91 1+, or = 9, (2.136) 

Net+Po,y = 91,2%1,2—-P1, yy I (2.13¢) 

No + Gort hy rr = 2(Pi,2t+ PF, y) — Pi ty (2.13d) 

3.1+P3,y = Pr.292,2+ (91,241 + 2,2) M,2 — P2,yy M1 — eaililia 7 ) 
.13¢ 


N3+ 931+ hg pr = $1,2P2,2+91,yP2,y+ (D1, cy P1,2 + P1, yyP1,y) I, 
+ hy e2(M1,2)*— P2,yt fect Pr, yy)” — Pr, yt Io (2.13f) 


At each step in approximation, the boundary conditions are now linear, all 
non-linearity involving lower-order functions which will have been determined as 
earlier approximations. All such terms have been placed on the right-hand side 
and are, at any given stage, known. 

As expected, the first approximation, ¢,, 7, are the solutions given earlier by 
Bellman & Pennington: 


?, = —/,sinh (4,t)e” cosx(= — By, ,e% cos x), (2.144) 
Ny, = cosh (2, t) cos x (= B,,(t) cosz), (2.146) 
where i= 1-2. (2.14¢) 


These solutions show the phenomena of cut-off wave-number since 4, = 0 at 
k = 1 and below this s, is positive so that surface waves grow. Note that for 
waves shorter than the cut-off wavelength (/ > 1), is imaginary and the waves 
oscillate; they are not damped. 

The second approximation 5, #2, extends the results published by Ingraham 
(1954) to include surface tension, although we solve only for a massless upper 


fluid: 


Byol(t 
d,= a sinh 24,¢+ ei sinh 2,t— st e?” cos 22, (2.14d) 
M2 = By,(t) cos 2a (2.14e) 
where ped = 2[1 — (2k)?], (2.14f) 
l (ta)? | fy)? | 
oo(t) = | sh 24,t- sh pot|). (2.14 
By, (t) 5 & a n= 1 Jeo h2u,t (7) cosh /igt (2.149) 
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The equation of the surface shows the growth of the first harmonic while the 
velocity potential shows the first harmonic and a purely time-dependent term 
which arises because of the additional pressure in the liquid required to produce 
the velocity of the rising fluid. The first harmonic is seen on analysis of these 
expressions to grow only by energy addition from the fundamental if the first 
harmonic wave-number is above the cut-off (i.e. the fundamental k > 4), while 
for first harmonic wave-numbers below the cut-off they grow also by direct 
energy addition from the acceleration field, once they are started. 

The third approximation introduces some new features. When the first and 
second approximations are introduced into the boundary conditions, equations 
(2.13e) and (2.13/), they become 


Ns. 1+ $3, y = Ps1(t) cos x + Py3(t) cos 32, (2.14h) 
N3+3,2+ kos, 27 = Q31(t) cos x + Q53(t) cos 32, (2.142) 


where the P and Q functions of ¢ are given in Appendix I. These equations show 
that as a result of the excitation of lower approximations not only a new mode 
appears (cos 3x) but also a feedback of excitation to the fundamental (cos x) 
occurs. The solution can be expressed as 


3 = C3,(t) cos x + C,,(t) cos 32, (2.15) 
Ns = B,,(t) cos x + B,,(t) cos 32, (2.16) 


where the B and C functions of ¢t are given in Appendix IT. The second harmonic 
behaves as the fundamental and first harmonic with respect to a cut-off wave- 
number above which it grows only at the expense of energy from the lower wave- 
number waves. The second harmonic is at the cut-off for the fundamental at 
} cut-off wave-number. Thus the distortion of the fundamental changes nature at 
each wave-number 1/n(n = 1,2,3,...) as each successive harmonic changes 
from direct growth to becoming oscillatory. 

The final solution for the growth of harmonic disturbances on the surface 
against air of an accelerated fluid is to the present order of approximation: 


n= ¥ 057, = No cosh “,t cos x + 92 Byo(t) cos 2x 
r=1 
+ 93(B,(t) cos x + By3(t) cos 3”) +.... (2.17) 
A corresponding sum holds for 
O(a, y,t) = y No Pr = — No ft, sinh p, te” cos x 
r=1 


t rk L F EB. 9 
+ 5 E sinh 2y1,¢ + a sinh 24,t- =.) e*” cos 22| 


+ 3(C3; cos x + Cy, cos 3x) +... (2.18) 
These series can be expressed as Fourier series by rearrangement: 
D 


g= > A 


—_ **m 
m=1 


(k, t) 46 and the B 


(k, t) cos ma, (2.19) 


where A,,(k,t) = > B -m(k; t) are given in Appendix IT up to 


— rm 
r=m 


B,,. Thus is displayed the harmonic content of the distorted growing disturbance. 


m 
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From the hyperbolic function nature of the B,,, functions already obtained 
and the manner of their formation, it appears that the entire series, if obtained, 
would change to a group of irrationally spaced time harmonic terms when the 
fundamental wave-number exceeds unity. This appears contrary to experi- 
mental results discussed in the next section. In any case, the present solution to 
the third approximation is still limited to small wave amplitudes; certainly not 
exceeding 7 = 7'k’ = 27, i.e. the wave amplitude equal to the wavelength. 
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Xesults of the analysis indicate that in order to get a reasonable representation 
of the wave-form even at comparatively early stages, at least a third approxima- 
tion is required. An example of the growth of a surface wave is given in figure 2 
where three consecutive shapes of the surface are shown for an initial sinusoidal 


disturbance of amplitude 9, = 0-3 and wave-number k = 0 (i.e. the effect of 


surface tension has been neglected). The asymmetrical growth of the surface is 
clearly indicated by the narrowing of the crest and the broadening of the trough 
of the wave with time. The influence of the amplitude and wave-number of the 
initial disturbance on its subsequent growth are shown in figures 3 and 4, 
respectively. As is to be observed there, asymmetry of the surface occurs much 
earlier for initial disturbances of larger amplitudes than for those of smaller 
amplitudes. Disturbances having wave-numbers close to ‘cut-off’ values (i.e. 
k = 1-0) are more stable than those having lower wave-numbers. Neglecting the 
effect of surface tension is equivalent to considering the limiting case of k = 0. 
Surface tension, therefore, has a definite stabilizing effect. Alternately, the 
limiting case k = 0 can be considered as a case of infinite acceleration for which all 


harmonics are below cut-off wave-number and are inherently unstable once they 
are started. 
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The present third-order theory also shows that, contrary to the prediction of 
a modified first-order theory (Bellman & Pennington 1954), the surface disturb- 
ance still grows at the ‘cut-off’ wave-number (i.e. = 1-0) in the manner given by 


N = No +5 Net? cos x + ge BI cos (2,/6t) — 1] cos 3x. (2.20) 


The crest and trough of the wave, however, are growing at equal rates. By 
further increasing the wave-number of the initial disturbance, as a consequence 
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FIGURE 4. Calculated wave profiles with different wave-numbers. 


k< 10 k= 5:0 k> 1-0 
Taylor’s original No exp (t) No exp (t) No exp (t) 
theory 
Bellman & Pen- No exp {/(1—k?) B No No COs {,/(k? —1) 
nington’s theory 
Third-order theory 7) exp {,/(1—k?) ¢ 3 mat? €5,(k) Noll (k? — 1) & sin {,/(k? — 1) * 


* Note. The term compared is the dominating one appearing in the amplitude A,(k, t) 
of the basic mode 4,(k, t) cos x. For e,,(k), see Appendix II. 


TABLE 1. Comparison between various theories regarding the growth of 
the basic mode in the earlier period 


of the feedback to the basic mode from the second harmonic, the surface disturb- 
ance, instead of growing exponentially in time as compared with those having 
wave-numbers k < 1-0, tends to be ‘overstabilized’, i.e. it oscillates in time with 
an ever-increasing amplitude. As a comparison between the predictions of the 
various theories, we have indicated qualitatively in table 1 the behaviour of the 
basic mode with variation of the wave-number. 
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B. Experimental apparatus 


The experimental results reported here were performed in an apparatus con- 
sisting of a glass-sided frame partially filled with liquid, a means of accelerating 
this frame and its contents, and various pieces of equipment for recording the 
behaviour of the liquid surface as the frame is accelerated. A schematic of the 
accelerator is given in figure 5. The frame which is accelerated is made of 
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FicurE 5. Arrangement of accelerator parts. 


aluminium with a glass front and back. Itsinsidedimensionsare 12in. x 5in. x lin. 
This frame is restrained to travel in a vertical direction by the two guide 
columns. 

At the start of an experimental run the frame is held at the top of the guide 
columns by a steel wire while tension is applied to the rubber tubes attached to 
the bottom of the frame. To begin the run an electric current is passed through 
the steel wire, melting it and releasing the frame. After the frame has left the 
area of interest another set of rubber tubes attached to the upper end of the frame 
comes into play and arrests the moving frame. 

In order to provide an initial disturbance to the surface of the liquid contained 
in the frame, a small paddle, its blade resting at the surface, is oscillated by a 
transducer driven by an audio oscillator and amplifier. 
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An indication of the initial amplitude of the travelling wave created by the 
paddle is obtained by observing the surface profile with a 30-power microscope. 
A stroboscope, flashing at a frequency close to the frequency at which the waves 
pass a fixed point, provides the illumination and creates the impression of a slowly 
moving wave profile which may be measured by means of a graduated reticle in 
the microscope. 

Photography is used for recording the time and position of the surface waves 
and also the position of the frame at millisecond intervals. A camera is placed in 
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FicuRE 7. Square root of frame displacement vs time. 

front of the frame, while behind the frame are positioned ten flash bulbs which are 
fired by switches tripped by the descending frame. Also located on the frame is 
a spark gap at which appears a spark discharge every millisecond. Because the 
camera shutter is open during the time that the frame is descending the film 
records: (a) the sparks giving the position of the frame every millisecond, (b) the 
ten positions of the interface at the times the flash bulbs were fired, and (c) the 
position of the spark gap at the time of the firing of the flash bulbs. The images of 
(b) and (c) above allow a precise determination of the relative times of the several 
events. A typical photograph is given in figure 6 (plate 2). 

The analysis of the film by travelling microscope and photographic enlarger 
provides the quantitative data. The travelling microscope is used to measure the 
position of the spark discharges, the wavelengths, and the wave amplitudes in the 
earliest stage of growth. The enlarger is used in the study of the later stages of 
growth. 

By plotting the square root of the distance from the initial position to the 
individual spark discharges against time (figure 7), the rate of acceleration and the 
relative time at which the acceleration began is obtained. 
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From extreme enlargements of the wave-forms of the later stages of growth an 
indication of the wave behaviour is obtained. A measurement of the distance 
from the wave spike or crest and of the bottom of the wave troughs from a mean 
horizontal line is made. The mean horizontal line is established by making the 
cross-sectional area occupied by the spikes equal to the area of the vacuated 


fluid in the troughs. This definition of a ‘mean horizontal line’ is used instead of 
a fixed line on the frame, because the meniscus growth changes the mean level of 


the waves relative to the frame. Figure 8 is an illustration of this later growth. 
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C. Experimental results on surface waves, small amplitudes 
The results of initial growth-rate experiments can best be discussed with the aid 
of figure 9 which is based on two non-dimensional numbers 7, a measure of the 
rate of wave-amplitude increase, as the ordinate and k, the ratio of wave-number 
to the cut-off wave-number, as abscissa. 

For k < 1 the predominate growth rate is that predicted by the linear theory of 
Bellman & Pennington (1954) and given here by equations 2.146 and 2.14c. This 
linear theory prediction is represented by a solid line forming a quarter circle in 
the wv, k-plane of figure 9. 

The value of from experimental results of wave growth in the k < 1 range is 
taken as the solution of the equation 


] 
fe = — cosh-! Im (2.21) 


> 
m No 


No is the ratio of the measured height to the initial wave height (for 


where 7),, 


experimental convenience the double amplitude or crest-to-trough distance is 
used) and ¢,,, is the non-dimensional time at which the wave reached the height y,,. 
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If the results of a large number of experiments performed by Lewis (1950) of a 
water-air interface are interpreted by the above method, then these results will 
fall in the shaded area of figure 9. As a result of the large accelerations and 
moderately long wavelengths used by Lewis the effects of surface tension would 
not be expected to be revealed. 

The single point placed by the names, Allred & Blount, represent the results 
of an experiment reported by them in which they used n-heptane and air as the 
two fluids. Unfortunately, this is the only experiment of this type reported by 
them; their other experiments involved the interface between two fluids with 
densities close to each other in value and, therefore, cannot be represented entirely 
as functions of the two dimensionless numbers used in this figure. 


Film k’ em-! g'/g k Nok, be 

D-3 8-64 2°32 1-00 0-11 0-49 
H-4 5-91 2-18 0-68 0-098 0-57 
M-7 6-64 2-24 0-76 0-11 0:64 
H-10 7-00 1-86 0-87 0-12 0-62 
H-11 7-26 1-96 0:88 0-024 0:46 
K-2 7°77 6-81 0-51 0-064 0-63 
L-1l 9-37 6-86 0-61 0-046 0-63 
L-3 6°53 6:94 0-42 0-032 0-61 
L-6 9-04 6°96 0-59 0-15 0-57 


TABLE 2. Initial wave-growth data 


The rest of the points in figure 9 are the results of the present work; the dimen- 
sional wave-numbers and accelerations of the experimental runs corresponding to 
the points in the #, k-plane are given in table 2. Although there is a definite 
decrease in growth rate as the wave-number approaches cut-off, all the photo- 
graphs of initially disturbed surfaces showed some evidence of wave growth. The 
experimental scatter is largely due to inaccuracies in the determination of the 
initial amplitudes. 

Beyond the cut-off wave-number as shown in figure 10 (plate 3) there is an 
indication of a type of behaviour here called over-stability. The photographstaken 
of wave growth under these conditions indicate a standing oscillation of the waves 
with amplitude increasing in time. In figure 10 the points A, B and C designate 
the same position in the space of the frame, and it is seen that the waves have 
been caught in different positions of oscillation. 

No more than a single oscillation has been observed; apparently after passing 
through one reversal of direction the wave grows to a large size and well-developed 
third-stage growth appears with the resulting competition between bubbles soon 
reducing the number of bubbles in a given space. There is some evidence that the 
larger the value of k for waves behaving in the over-stable manner, the higher the 
frequency at which the wave oscillates, as would be expected from equation 
(2.14c). 

D. Later stage wave-growth experiments 
The films which were analysed for third-stage growth showed considerable 
scatter of results even though a selective process was applied in the choice of 
which group of waves were to be studied; i.e. those waves were chosen which 
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appeared to be in a group of waves of a fairly uniform nature. The scatter was 
largely due to two things: the presence of odd modes of vibration imposed on the 
steady growth, and to small variations in initial conditions which become 
magnified at later times. 

“xcept for the prediction of over-stability, the third-order theory does not 
check well with the experimental results. The theory predicts an ever-growing 
oscillation while the experiments never show more than one cycle of oscillation 
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before monotonic growth occurs. It is not known at present whether the in- 
adequacy is in the theoretical treatment or in the experimental boundary 
conditions. 

Spike acceleration is illustrated by figure 11. Of the photographs measured for 
spike growth all indicated an acceleration of less than half the rate of free fall in 
the particular virtual gravity field applied. This may be compared with Carrier’s 
(1953) work and the machine calculations of Pennington et al. (1953), which 
indicated a rate of acceleration equal to the free-fall rate. Since both of these 
theoretical predictions assume the surface tension to be zero, they would indicate 
a higher rate of growth because the surface tension forces act against the direction 
of acceleration of the spikes. 

On the other hand, experimental results on bubble velocity resulted in good 
agreement with the theoretical predictions. The results, which are illustrated by 
two examples given in figure 12, in general give the conclusion that the value of 


c in the equation, v = c, gA, lies between 0-2 and 0-3. This agrees with the theory 
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of Birkhoff (1954) and the work of Garabedian (1957), even though they did not 
consider surface tension effects in their analyses. 

Examples of competition between bubbles are observed in every photograph 
taken of the wave disturbance type. 
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FigurRE 12. Bubble velocity. 


3. Experiments involving a liquid sheet at the wall 

If the surface is not disturbed by the vibrating paddle, a distortion of the 
surface still results from the liquid meniscus around the edge of the accelerated 
container as is illustrated by figure 13 (plate 4). Point A is the position reached 
by a very thin sheet of liquid next to the glass. Point B is the top of a thin wall of 
liquid that is in the exact centre between the glass walls. Points C and D are the 
tips of two two-dimensional bubbles which occupy the space between the sheets 
at the wall and the liquid wall in the centre. The more or less round bubbles of 
which there are four, one in each corner, are at the points E and F. The above is 
further illustrated by figure 14 (plate 5) in which an oblique view of the same 
type of phenomenon is offered. 

Some study was made on the behaviour of the liquid sheet at the side wall. 
The results of this study are given in figure 15. Two liquids, methanol and carbon 
tetrachloride, were used in these experiments, and little difference was found 
between the two except that in the case of the carbon tetrachloride a second 
thickening of the sheet was seen to occur some distance back from the leading 
edge. This thickening was seen to catch up with the leading edge, after the leading- 
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edge acceleration began to slacken, and to provide fresh impetus allowing accelera- 
tion to continue at the original rate for a further distance. This point is illustrated 
by a comparison of figures 15 and 16. In figure 16 the line, labelled second wave 
line, represents the region where the thickening has caught up with the leading 
edge. 

A graph, figure 17, of the ratio of the leading-edge acceleration to the applied 
gravity field versus the applied gravity field indicates that the leading edge 
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FicurE 17. The ratio of the acceleration of the leading edge of the liquid sheet at the wall 
to the acceleration of free fall in the applied gravity field vs the applied gravity field. 














Taylor instability of finite surface waves 








0} A 90 A of 7) A 
3 7) Ae 
= | i |t | lt 1} 
5 
52 | 2} | 2 2} 
= 3 3} | 3] 3+ CS 
8 } | | | 
= 4 ;4¢ | 44+ | 4 
™ ; | | | | 
(a) ere Tee } 51 5 2 
2 l | ae r @2 FF (G2 | ae | 
Time (msec) 43-7 60-0 72°8 83-3 
Displacement (em) 1-17 2-22 3-30 4-30 


FicuRE 18a. Liquid sheet profiles at a low acceleration rate. 
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Figure 186. Liquid sheet profiles at a high acceleration rate. 
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FicurE 19. Profile thickness. 
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travels, in the ranges of applied gravities used here, at from 0-6 to 0-7 as fast as 
a free falling body. This is a higher acceleration than any observed in the case of 
wave spikes. 

In order to get more information on the manner in which the accelerating sheet 
acts, the profiles of the sheet at the wall were determined at the highest and 
lowest rates at which the experiments with methanol were run. These results are 
given in figures 18a and 1856. 

From the profile, a characteristic thickness can be measured and plotted versus 
time as in figure 19. The points marked A and B in this figure represent the time at 
which the displacement of the leading edge no longer maintained its constant 
acceleration. From this it can be deduced that at the time A for the experiment 
giving the profile shown in figure 18a and at the time B for the experiment giving 
the profile shown in figure 18 the liquid sheet ran out of fluid in the sense that the 
rate of supply from the main body of liquid was insufficient to allow things to 
continue as they had been. 


The authors gratefully acknowledge the support of this work by the Office of 
Ordnance Research, U.S. Army, through Combustion Research Contract 
No. DA-19-920-ORD-1029. 


Appendix I 

Amplitudes of forcing functions of the third approximation, equation (2.14h, 7) 
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FicurE |. A later stage of wave development. 
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FicureE 6, Typical photograph of position of interface at ten stages of growth. 
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Figure 10. Example of over-stable growth. 
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Plate 3 





Journal of Fluid Mechanics, Vol. 7, part 2 


—_—— 
a 
meee 


op hety! 
“4 " i 

| r i Y 
lpi 


at tener, 
" 


Freureé 13. The growth of an initially undisturbed surface. 
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Ficure 14. Oblique view of the growth of an initially undisturbed surface. 
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Appendix II 
Amplitude functions of the third-stage approximation, equations (2.15) and (2.16) 
Cg, (t) = Poi(t)— Bai(t), Cga(t) = Pya(t) — 3 Bga(t), 
B5,(t) = ag, cosh (4, t) + bg, cosh (34, t) + cs, cosh (4, ¢) cosh (6) 


+d3,(7) sinh (4, ¢) sinh (gt) + es) 4, ¢ sinh (j2,¢), 
B33(t) = agg cosh (/14t) + bg cosh (3/4, ¢) + ¢33 cosh (4, t) cosh (st) 
+ da “1) sinh (j;¢) sinh (j1,t) + eg5 cosh (Hs), 
He 


where P,,(t) and P,,(t) are given in Appendix I, 2 = r(1—r?k?), and 


1/3 1-yj 5 (4/9)? — 1) 
ee ed |) ee b. _ ot Pll’? 
a31 (63: +€31), 93, 8 | 39 _ mn ili) le a ATE 
—— ails :) 2(/4,/ Hs) 
= ds, = 
~ [4(Hy/42)? — 1)?’ a (4/2) se 


9 1—pj l ni Hs Pett )?— i. 
Cy = ’ 


64 16 fal Me oa, l 


l (271-,72 4( 0, /f5)? + 1) 
a3g = — Bi ea > , 
1 — (fl3/,)* (32? na 
ee ae (9 1=—KG 3 4(ay/My)? +1) 
at 9 — (fg/fty)? (32 8 eae /y/ to)? — VY? 
a" (44/2)? ] r 
- © 4( tty / fly)” — 1 (Mg/fey)* — 2 [1 + (Meo/ ea)? (Hg/ a)? + = (Mee)? ? 
dy, = —6 — Hala)” 1+ (H2/y)? — (Hs/#a)” ren 
Z 4 (fy) fa)? — 1 (Mg/y)* — 2[1 + (e/a)? ) (Ms/a)? + L — (Hel) P 
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Dispersion due to molecular diffusion and macroscopic 
mixing in flow through a network of capillaries 


By P. G. SAFFMAN 


Cavendish Laboratory, Cambridge 
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This paper is concerned with the dispersion of a material quantity in the steady 
flow of a viscous fluid through a random network of capillaries (which is a useful 
model of a porous medium), for the case in which molecular diffusion and macro- 
scopic mixing, due to the randomness of the streamlines, are both important. 
A Lagrangian correlation function is introduced and the longitudinal and lateral 
effective diffusivities are thereby calculated for all values of Ul/« less than some 
large value. Here, / denotes the length of a capillary, U the mean velocity of the 
fluid, and x the molecular diffusivity of the material quantity. The theory is 
compared with experimental observations of dispersion in flow through granular 


beds. 


1. Introduction 

In a recent paper (hereafter referred to as I), Saffman (1959) has investigated 
the dispersion of a dynamically neutral material quantity in a viscous fluid 
flowing through a porous medium, on the basis that the porous medium may be 
regarded as equivalent to a statistically isotropic network of straight capillaries. 
The discussion was restricted to the case in which the dispersion is primarily due 
to macroscopic mixing. This arises from the randomness of the streamlines 
through the network, and has an effective diffusivity of order Ul, where I is the 
length scale of the capillaries and U is the mean velocity of the fluid (the inter- 
stitial velocity). The results of I are then only valid when 


k < Ul, (1.1) 


where « is the molecular diffusivity of the material quantity. 

The method of I utilized the fact that the path of a fluid particle through 
a random network of capillaries can be treated as a random walk of randomly 
orientated steps, the velocity in a step depending upon its direction. This approach 
is clearly valid only if (1.1) is satisfied. An alternative approach, involving the 
idea of Lagrangian correlation functions which are well known in the theory of 
turbulent diffusion (see § 3 below), was suggested to the author by Dr T. H. Ellison 
in a private communication. As it happens, this method seems to be in general 
(although not always?) less useful than the random walk approach for the cases 

+ For example, it gives the dispersion in certain cases for small values of the time from 


the initial instant, whereas the random walk approach is an asymptotic theory for large 
values of this time. 
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considered in I. However, it affords a means of calculating the dispersion when 
(1.1) is not satisfied, i.e. when the effect of molecular diffusion is comparable with 
or greater than the macroscopic mixing, and this is the purpose of the present 
paper. The theory is not without practical importance since a network of capil- 
laries is not an altogether unreasonable model of a porous medium and actual 
flows in which Ul/« < O(1) are not uncommon, e.g. the flow of gases through 
porous media. 


2. The flow in the network 

The work is restricted to the case in which the mean (macroscopic) velocity in 
the network is uniform and unidirectional, and the Reynolds number of the flow 
in the capillaries is sufficiently small for viscous effects to be dominant (this is 
equivalent to supposing that the mean flow satisfies Darcy’s law). We also 
suppose (primarily for simplicity, as in I) that the pores are circular of radius a 
and of equal length 1, where a </. Then to a reasonable approximation, the 
average velocity in a capillary making an angle 0 (0 < 0 < 47) with the direction 
of mean motion may be taken as 


q = 3U cos8, (2.1) 


where U is the average velocity of the entire fluid (see I, equation (2.19), where 
q is denoted by 4). 
Since this paper is primarily concerned with the dispersion when k is comparable 

with or greater than Ul, we shall suppose further that 

Ul _ SF i. @ 

x " a @ > 8k" = 
Now a?/8x is the time scale for molecular diffusion to smooth out variations in 
concentration across the cross-section of a capillary, and the condition (2.2) 
expresses the condition that this is small compared with the time spent by a fluid 
particle in a capillary. The work of Taylor (1953; 1954) and Aris (1956) shows that 
in this case the concentration is approximately uniform over the cross-section and 
disperses, relative to axes moving with the velocity q, according to a diffusivity 


a*q? 
K+—— =D, say. 2:3 
a (2.3) 


Thus, when (2.2) is satisfied, the dispersion in the network is the same as if the 
velocity in each capillary were uniform throughout it, with the value q, and the 
material is subject to a diffusion process along the capillary with a diffusivity D 
(which depends upon the direction of the capillary). Note that 


aq? 3 a2 =) i 3 Ul 


48x? > 162? 74 


‘ey  2e 
by virtue of (2.2); so that D = x if Ul/« is not large compared with unity, and for 


all Ul/k if a €l. 
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3. The Lagrangian correlation function 

If u(t) is the velocity component in the z-direction of a marked particle 
participating in a random process, the Lagrangian correlation function or 
covariance of the velocity u at time ¢ and at a later time ¢’ is u(t) u(t’), where the 
bar denotes an average taken in an appropriate way (which will depend upon the 
precise nature of the random process). As was first pointed out by Taylor (1921), 


*t 
the displacement in the x-direction after time ¢ is X(t) = | u(t’)dt’, and 
dX /dt = u(t). Hence shes 
l= ; ‘ ; 
5 X2(t) = 2X (Hult) =2 | ult)ut’)dt’. (3.1) 
¢ J0 


For the case in which u(t) is a stationary random function of the time 
(e.g. turbulent diffusion along a pipe; see Batchelor & Townsend 1956), 
u(t) u(t’) = R(t—t’), where F(z) is the correlation function of the velocity w for 
a time interval 7 (and depends only on 7 and not on the time from any given 
instant), and the axes are chosen so that « = X = 0. The overbar here denotes an 


ensemble or stochastic average. Then (3.1) becomes 


dX? ‘ys 
— D 99 
as 2) R(r) dr (3.2) 
°t 
and hence X?=2 | (t—7) R(r) dr. (3.3) 
/0 


Thus when ¢ is large compared with the time scale for which the velocities at 
successive times are effectively correlated, 


here a. 


X2 x2 | R(r)dr—2 | rR(r)dr, (3.4) 


0 70 


assuming that the integrals converge. 


X2 Proc 
Thus, ho | R(r) dr, 
~ Jo 
Par 
and hence | R(r) dr may beregardedas the effective diffusivity. Itshould perhaps 


/0 
be pointed out for the benefit of readers unfamiliar with this analysis that (3.4) 


does not prove that the dispersion of a cloud of marked particles is given (asymp- 
totically) by solutions of the diffusion equation with this effective diffusivity. To 
do this, it is necessary to prove that the probability distribution of X is normal, 
and that the displacement of different fluid particles in the same realization of the 
flow become statistically uncorrelated. However, there are usually no reasons for 
believing that these requirements are not satisfied, even though in most cases the 
formulation of rigorous proofs appears to present great difficulties. 

If now u(t) denotes the longitudinal component (i.e. in the direction of mean 
flow) of the velocity ofa marked particle in steady flow through a random network 
of capillaries, we may proceed in an identical manner, provided the network is 








be 
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statistically homogeneous. For convenience, we refer uw and X to axes moving 
with the mean velocity, so that 


u = (q+ Qp)cosd—U = U(3 cos? 6—1)+qpcos8, (3.5) 


where 0 = ((t) refers to the capillary containing the particle at time ¢. qp is the 
random part of the velocity along the capillary due to the diffusion process with 
diffusivity D; it can be regarded as the velocity due to a Brownian motion 
equivalent to the diffusion process. That is, the velocity of a marked particle 
along a capillary fluctuates randomly because of the actual Brownian motion 
(caused by the collisions with other molecules) and because, as it moves to and fro 
across the capillary cross-section, the local mean velocity varies. These effects are 
mathematically equivalent to diffusion with diffusivity D and a corresponding 
fictitious Brownian velocity ¢, with the properties 


. 


(Ip) =9 and | <qp(0) qp(7)> ar = D, (3.6) 
0 


where ¢ ) denotes an average with respect to the Brownian motion. 

This argument is, of course, valid only if the time interval over which successive 
values of gp are correlated is small compared with the average time spent by 
a marked particle in a capillary. The former is at most a?/8« (and will be much less 
if molecular transport dominates the diffusion process), whereas the order of 
magnitude of the latter is the lesser of 1/U and /?/2«, and the argument is therefore 

valid by virtue of (2.2) and the condition a < l. 

For the case of turbulent flow, the average denoted by the overbar is an ensemble 
average over all realizations of the flow. For the case of steady flow in a random 
network, the average is over all configurations of the network, but this is equiva- 
lent, owing to the statistical homogeneity of the network, to an average over the 
initial position and velocity of the particle. Thus, if f denotes a property of 
a marked particle which at timet = 0wasina pore making an angle @(0 < 6 < 47) 
with the direction of mean motion and with azimuthal angle ¢ (0 < ¢ < 27) about 
the same direction, was at a distance 1/5 (0 < € < 1) from the entrance of the 
capillary, and had a‘ Brownian velocity’ 7p, where the values of £, 6 and #7 = cos@ 
are uniformly distributed (since the proportion of pores in the range @ to 0 + d@ is 
sin #d@), then d 


» 


Pon ff 


fl (2 1 
| fsin 0 dé dé 
“0 /020 I 


inf” 


a) is 
7 oP (dp): (3.7) 








where p(qp) is the probability that the ‘Brownian velocity’ is less than gp). Note 

that sie 

f) = | fap (qd). 
/0 


The value of f is independent of the precise origin of the time. Moreover, 
because of the statistical homogeneity, u(t) = u(0) = 0, and u(t) u(t’) is afunction 
of (t’—t) only, equal to R(t’ —t), say. Thus the covariance is given by 


R(r) = u(0) u(r), 


with the dispersion following from (3.3) and (3.4). 
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4. Calculation of the covariance 


Let primed quantities denote values after a time 7 from some initial instant. We 
have, writing “7 = cos@, yz’ = cos6’, 


uu’ = U?(3u? — 1) (342-1) +U'InIp t+ Ufugn(3e'2— 1) + 4'qp(3u?—-1)}. (4.1) 


The last term in (4.1) is zero since (q¢p) = 0 and » is independent of gp. 
Consider now the second term. The time interval over which qp and qj are 
effectively correlated is small compared with the average time for the particle to 
leave the initial capillary and enter another. In other words (qgpqp) is non- 
negligible only for values of 7 so small that the value of « does not change (except 
for a negligible range of values of for which a particle is next to a pore junction). 
Hence, to a reasonable approximation, 
, aPamtees Ponen erenne (C1 
MEDI = 4p dp = | GDI) de 
on integrating with respect to £ and ¢. (This is permissible since (qpqp) will 
depend, to the present approximation, only on 7 and the mechanics of the 
diffusion process with diffusivity D.) 
It now follows from (3.6) that 


ee 1 
[ MK'GpIpdt = | DrPdu (if t > a?/8x) 
Jo Jo 
Kk 3 aU? 
— _ ara Se 9 
a 3° 80 K a 


on substituting from (2.3) and (2.1). 
Thus the problem reduces to the calculation of the first term in (4.1), namely, 


(3u2—1)(3u’2—1) = S(7), say. 


There are two distinct contributions to S(7). First, there is the integral over 
those initial conditions such that the marked particle has not reached either of the 
ends of the initial capillary after time 7, so that ~ = y’; this gives a contribution 
S,(7), say. Secondly, there is the integral over the conditions for which the particle 
is in a different capillary after time 7 (or has returned to the initial capillary after 
having passed through one of the junctions at its ends), giving a contribution 
S,(7), say. Then S(7) = S,(7) + S,(7). 


4.1. The calculation of S,(7) 


Let P(&9.,7) denote the probability that a marked particle released at time 
T = Oata distance /f, from the entrance of the capillary has not reached one of the 
ends after time 7. Since « = yw’ for a particle which remains in the initial capillary, 


' ine dd 
S,(7) = | (34° — 1)? Pdudgy5—. (4.3) 
JoJoJo on 





(P(&o. 4,7) is in fact the integral of dp(q,p) over those modes of the ‘ Brownian 
motion’ such that £ + 0 or 1 in time 7.) 
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To calculate P, we note that the marked particles are convected along the 
capillary with velocity q (= 3U,) and diffuse relative to an origin moving at this 
speed with a diffusivity D. If they reach the ends of the capillary, they no longer 
contribute to S,(7). Hence the probability that the particle lies between / and 
l(€ + dé) after time 7 (not having passed through £ = 0 or 1) isc dé, where c(&,7) is 
the solution of ‘ . 
ce qcee D 


¢ 
eae re (4.4 
or Lae ae? ie 
subject to the initial condition ¢ = 6(§—£,) when 7 = 0 (4(&) is the Dirac delta 
function) and the boundary conditions c = 0 when € = 0 and 1 for all 7. Then 


1 
P(£o, 4,7) = [ cdé. (4.5) 


J0 
The solution of (4.4) is straightforward. It is easily verified that 


io 8) 
C = Ya, sin nngé eS exp { — (n?22 + M?) Dr/l?}, 
1 
l 2Ulu 
where Buen a, 4.6 
; 2D K+3a?U2u?/16K- 9 


is a solution of the differential equation which satisfies the boundary conditions. 


It satisfies the initial conditions if 


@ 
6(§ —£,) = Na, sin nag e™, 


from which it follows by the usual methods of Fourier series that 


a, = 2e-“osin n7€,. 


n 


Then (4.5) gives 


fo) 9 
2n7z 
P(&q, ¢,T) = ——— 
So» / 7 n2n? + M? 


4 


—(—1)"e} eosin n7é, exp { — (n?2? + M?) Dr/l*}. 


(4.7) 


f 
L 


On substituting in (4.3) and integrating with respect to £) and ¢, we obtain 


Sr) =f 4du(32—1yes 
en * BN 7 (n?7? + M?)? 


x {1—(—1)" cosh M} exp { —(n?7?+ M?) Dr/l?}, (4.8) 


where, in general, D and M are functions of ~ = cos@. 
This expression is fairly complicated. However, for the purpose of finding the 


effective diffusivity we only require the value of [ “S,(r) dr. Integrating (4.8) term 
/0 
by term and summing the infinite series, we find 
oo 1 co f2 n272 
S,(7)d7 = [ 4du(3n? — 1)? > — {1 -(—1)" cosh M} 
[ 1( ) Te j ( I ) 7 D (n2n? + M?)3' ( ) J 


' l Ea 2 
- ae ee besa 
-{ dt op BH 1) (;,coth M in). (4.9) 


0 
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We have here made use of the following formulae (which may be obtained by 
expressing cosh Ax as a Fourier cosine series in 2 and putting x = 0 and 7): 


coshA—1 | holt Saad Toe) cosh aA 


nn? +22 


A—sinh a _ A al 


242sinhA = T n2n2 + A?’ 


and the formulae which can be obtained by differentiating them with respect to A. 


4.2. The calculation of S,(r7) 


The determination of S,(7) is rather more difficult. We shall assume, first, that 
the values of ~ for successive capillaries are statistically independent; that is, 
when a particle reaches a junction the choice of the next capillary is independent 
of the previous one. The larger «/Ul the better will be this assumption, because 
molecular diffusion has a much smaller correlation time scale than the macro- 
scopic mixing. 

It might now be thought that 


However, this argument is false because it ignores the fact that although y’ is 
independent of , it is not independent of the fact that the particle has passed 
through a junction during the time 7 and this, as will be seen below, affects the 


(°20 


probability distribution of ’. Nevertheless, it is in fact true that | S,(7)d7 = 0, 





J0 


he.o) 
so that S,(7) does not contribute to the effective diffusivity. | 7S,(7) dr + 0, 
70 
so there is a contribution to X*. The remainder of this section is devoted to the 
calculation of S,(7), or, to be more precise, its Laplace transform which is 
sufficient for our purposes. 

Let Q(u,t) du denote the probability that a marked particle released at a junc- 
tion at ¢ = Ois in a capillary in the range uw to ~+dy after time ¢. 

The probability that a particle released at t = 0 at a distance lg from the 
entrance to a capillary with a given value of reaches a junction in the time 
interval from ¢ to t+dt is ‘ 

= {1— P(E, nu, t)} dt, 
where P(£,,f) is given by (4.7). Hence the probability that this particle is in 
a capillary in the range yw’ to ”’ + dy’ after time 7 from the initial instant of release 
is ‘ 
dy’ | O('.7 —t) = {l—P(é, 4, t)}dt = L(E,4.4',7) du’, say. (4.10) 


70 


Then by definition 


l/l 
S,(7) = | | dé du dy’ (3p? — 1) (3? — 1) L(E, wy’. 7). (4.11) 
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If all the pores leading from a junction had equal probability of being entered 
by a marked particle, then Q(u, t) = 1 and S,(r)= 0. However, this is not the case 
since there is a tendency for particles to enter pores in which q is large rather than 
those in which it issmall. (If« < Ul, the probability ofa particle leaving a junction 
by a pore in the range of width dy is 2u dy, and Q(u,t) = 2 fort <1U.) On the 
other hand, for large values of t the particle will (so to speak) forget where it was 
at ¢ = 0 and the probability of the particle being in a capillary in a particular 
direction will become independent of the direction, i.e. Q(u,t) > last —> oo. But 
since S,(7) depends upon the values of Q(u,t) for all values of t between 0 and 7, 
S,(7) will not in general vanish. 

It is useful now to take Laplace transforms with respect to the time. These will 
be denoted by an asterisk, i.e. the transform of f(t) is f*(p). Then since (4.10) is 
a convolution, we have on integrating with respect to €, 


rl 
2(p) = < [ dy du’ (3n? — 1) (342-1) Q*(u’, p) {1—pPF(u.p)}, (4-12) 


~0 0 
4n?7* 1—(—1)"cosh M 
where Pr(“#,p) = > , 4.13 
iP) T (n?27? + M?)? p+ (n?20? + M?) D/P 
Now for the effective diffusivity, we really require only S,(7) dr and this is 
J0 
S}(O). It now follows from (4.11) and (4.12) that S$(O) = 0, provided 
“1 
lim | (342-1) Q*(n’, p) dy’ (4.14) 
p—>0 J/0 


is finite. We know intuitively that Q(«’,t) > 1 as t > «, and if we assume that 
Q(y’, t) — Lisintegrable from 0 to «, it follows that Q*(u’, p) = p-'+ O(1) asp — 0, 
and hence that (4.14) is finite. However, this last assumption is far from obvious 
(Q(u’,t)-—1 might be O(1/t)) and it therefore seems worth while presenting an 
explicit calculation of Q*(y’, p). 

Let E(y’,t)dy’ denote the probability that a particle released at a junction is 
after time ¢ in a capillary leading from this junction in the range yp’ to yz’ + dy’, and 
in this time ¢ has not returned to the junction or gone to the other end of the pore. 
The value of y’ in (and only in) E(y’,t) ranges from —1 to +1 (elsewhere 
0 < yw’ < 1),totake account of the fact that the particle may diffuse into a capillary 
against the direction of the stream. Then Q(w’,t) is related to E(w’,t) by 


hd 


rt “ 
O(m’.t) = E(w’,t)+ B(-w’,t)+ au Qu’), |= | But yan, (4.15) 
0 ¢ 


e 


-1 


the last term being the probability that it was at a junction sometime during the 
interval! 0 to ¢ and went on from there into a pore with ’ in the appropriate range. 
Taking Laplace transforms of (4.15), we obtain 


ee | 
Q*(u',p) = B*(u', p)+E*(—w'. p)+Q*(u’ pt —p | Bu. pda, (4.16) 
1 
P+] 


since E(u, 0)du = 1 


7-1 
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7 E* a . E* = Zs 
Hence Q*(n’,p) = a i—,% p) (4.17) 
p E*(u,p)du 
ot 


Now E(u, t) may be determined in a manner similar to that for P(£, , t), except 
that care is required because the particle starts at a junction. We overcome this 
difficulty by taking c = 1/A7’, for € = 0 and 0 < t < AT’, as the boundary condi- 
tion on the concentration c(£, t) along the capillary, where c(é, t) satisfies (4.4); and 
c = Oinitially,c = Oat £ = 1 forallt,andc = Oat = Ofort > AT. Then, we take 

a - 


E(u,t) = lim c(f,t)dg, and E*(u,p) = lim c*(£, p) dé. 
AT->0 J0 AT—0 J0 


Taking the Laplace transform of (4.4), we find that c*(&, ) satisfies 


o2c* oc* [Fp . 
—2M ———c* =0, 18 
ag f 36 D Cc ) (4.18) 


which has to be solved subject to the boundary conditions 
] " 
c*(1,p)=0, c*(0,p) = part —ePAT} +1 as AT—0. 
The solution of (4.18) is perfectly straightforward (it is a linear equation with 
constant coefficients) and we find after a little algebra that 
M sinh (M? +[?p/D)* + (M? + 2p/D)3 {cosh (M? + I2p/D)t — eM} 


E*(,p) = 
(4, P) (I2p/D) sinh (M? + [2p/D)} 


(4.19) 


(4.19), (4.17), (4.13) and (4.12) thus give an explicit expression for S}(p), and 
hence determine, in principle, S,(7). We find from (4.19) and (4.17) that, as p -> 0, 


,~.[l|-MecothM ff}, 1—McothM ; 
pO*(u,p) = 1+ ip 12D a du wp | (4.20) 
thus verifying that (4.14) is finite and that 
S,(r)dr = 0. (4.21) 


#0 


5. The longitudinal dispersion 

Collecting together the results of §4, we have from (4.2), (4.9) and (4.21), the 
result that the longitudinal dispersion in a network of capillaries is asymptotically 
described by an effective diffusivity 


x 2240 FOr , Uf coth M—1 





kK, = ea ae - “ef (3° — 1) Da? du. (5.1) 

3 Ulu 302 U2 2 
7} are J ~—-4 — :: ) =K wy ai 5.2 
where I 3D and D=K+ —* (5.2) 


The expression for x; is valid provided we may replace the upper limit of inte- 
gration in the expression (3.3) for X? by t = 0. This is permissible when the 
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exponential term in (4.8) involving the time is small compared with unity and 
S¥(1/t) = S$(0), which is the case if t > 1?/x. 

When there is no convection and no macroscopic mixing (i.e. U = 0), it follows 
from (5.1) that x, = 4«. Now experiments on dispersion in granular beds (e.g. 
Carberry & Bretton 1958) give values of x, = %x, when convection is negligible. 
The discrepancy seems to be due quite simply to a granular bed not being 
a random network of capillaries. 

However, this is not to say that the present theory is not applicable to disper- 
sion in granular beds; it is at least of qualitative value. Moreover, when there is 
flow through a granular bed, the streamlines will approximate to those through 
a random network of capillaries, and it is not too unreasonable to suppose that 
the macroscopic mixing due to the meanderings of the streamlines is given 
approximately by the terms in (5.1) involving U. In other words, we suggest 
that (for the purposes of comparison with experiments in granular beds) the 
above expression for x, be replaced by 


2772 1 M-— 
3 aU | M coth M hay, (5.3) 


Kr = Kn +35, + 112U? ; (3? — 1)? — DI 
where x,, is the measured diffusivity when U = 0 (and M and D are given by 
(5.2)). The closer the geometry of the bed approaches that of a random network, 
then presumably the closer will be x,, to 4«. The calculating of x,, presents great 
difficulties; however, it may be readily measured by means of an electrical 
analogue (see, for example, Wooding 1959). 

The expression (5.3) for x, still contains two quantities, a and 1, whose values 
have to be fixed in some way from the properties of the granular bed. For a, 
areasonable value isa = (24k/c)!, where k and o are the permeability and porosity, 
respectively, of the bed. The bed and the random network will then have the same 
permeability (see I, §2). In any case, the values of x, are fairly insensitive to 
changes in a, provided a/l is not greater than about } (which is about the largest 
permissible value consistent with the assumptions about the flow through the 
capillaries). 

As regards 1, probably the best value to take is the mean diameter of the grains, 
and it is not possible to be more precise than this. However, it may be an under- 
estimate, especially if there are cracks or fissures or channels running through the 
bed. The presence of these would be equivalent to increasing the value of 1, since 
they roughly correspond to longer capillaries. We can in fact work out the theory 
for a network of capillaries of varying lengths; this simply gives (5.1) averaged 
over the lengths of the capillaries, but since the distribution of capillary lengths is 
unknown, there is no real loss in supposing them to be all equal to some suitably 
chosen value which may be somewhat larger than the mean particle diameter if 
channels run through the bed or the packing is not too uniform. 

The longitudinal diffusivity as given by (5.3) may be evaluated numerically for 
given values of a/l and Ul/«. For Ul/k < 1, 

1 U® 


5.4 
m 15 K (5 ) 
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and for Ul/k > 1 (but not so large that (2.2) is violated), it can be shown that 


_ Ulf SOL AT la? Ul) KK?) ay) 
kK, = 6 jog. “To-8F x +k +$e+0(75] ; (5.5) 


The expression (5.5) may be compared with the equivalent expression obtained 
in I (equations (4.42) and (4.5)) by means of the random walk analysis, namely 
Ul 3Ul 1 


ar (5.6 
i“ 6 2K 6C 


log, 
The leading terms agree, providing a useful check on the analysis. The difference 
in the other terms seems to arise from the fact that in I molecular diffusion was 
taken into account only in an approximate manner. 
Hiby (1959) gives the following empirical formula for longitudinal dispersion in 
beds of glass spheres: ee 
; ae 0:65 Ul oa 
K, = 0-67K 4+ =e (0.7) 
1+6:7(k/U1)2 
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Figure 1. Comparison of theory and experiment for the longitudinal dispersion. 


, Hiby (1959); - - - - - . Blackwell et al. (1959); x, Carberry & Bretton (1958); 
, calculated values from (5.3) with a/l =} and k,, = 3k. 


where / is the mean diameter of the spheres. According to this expression, 
macroscopic mixing is completely negligible for Ul/K < 0-5. On the other hand. 
Blackwell, Rayne & Terry (1959) suggest on the basis of experiments in beds of 
sand the empirical expression: 
kK; UN ; ; 
= s'8| for > 0-5, (0.8) 
K K 


K 


where / is the mean particle diameter. Also, they find that macroscopic mixing 
becomes significant at about Ul/x = 0-04. 
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These results are shown on figure 1, together with the theoretical curve 
taking a/l =} and x,, = §x (a/l = 75 is probably a better value for a bed 
of sand but the difference in the calculated values is slight). It will be seen 
that the theory predicts that macroscopic mixing is negligible for Ul/k < 0-5, 
and the calculated values are in remarkably good agreement with the expression 
of Hiby. Nevertheless, the discrepancy between the experimental observations 
is such that it is not possible to say definitely whether or not the theory is 
valid when applied to granular beds. A possible explanation of the differences 
might be that the packing of the sand in some of the experiments was not 
uniform, in particular there may have been gaps or channels along the walls of 
the cylinders containing the medium. Hiby states that he was unable to obtain 
uniform packing right up to the walls of his apparatus, and he went to great 
trouble to eliminate the wall effect which could affect the results to a factor of two 
or more. As mentioned previously, non-uniform packing or the presence of 
channels along the walls would be equivalent to an increase in the value of 1, and 
would thereby reduce the discrepancy between the calculated values and the 
empirical curve of Blackwell et al.; but it is not yet possible to say whether this 
explanation is sufficient. 


6. The skewness of the longitudinal dispersion 


In the previous section, it was tacitly assumed that the longitudinal distribu- 
tion of concentration was Gaussian, so that it is meaningful to regard X?/2t as an 
effective diffusivity. We shall now obtain an estimate of X% in order to see how 
Gaussian the distribution actually is. This calculation is not without practical 
importance since observed concentration profiles are often skew. 

By arguments similar to those in § 3, it follows that 


d ee ce 
X3 = 3X%(t)u(t) =3 | | wt) a(t?) w(t") ae’ ae” 








it 0/0 
re “ws 
=@6] dr, R,(71, 72) 47; (6.1) 
/0 70 
where R,(71. 72) = u(O) u(71) u(T2) (6.2) 


is a triple Lagrangian correlation function. 

Now there are two types of contribution to #,(7,,7,). First, there are the con- 
tributions from the cases in which the particle is in the same capillary after time 
7, from the instant of release; and secondly, from the cases in which it is in 
a different capillary after time 7,. We shall assume that the second contribution 
has negligible effect on the value of X4(¢) for large values of t (as is the case for X?). 
Although this assumption is highly plausible, a rigorous proof seems to present 
great difficulties. 

To calculate the first contribution, we require the probability that a particle 
released at a given position is in the same capillary after time intervals 7, and 7,, 
where 7, < 72, and the particle has not been to one of the ends of the capillary and 
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come back again (this is counted as being in a different pore). This probability is 
just P(g), 4,72), a8 given by (4.7). Hence, for 7, < 75, 
el fl 


Rylryt2) = U8 | | (32-1) Pop, 72) dE 
0/0 


Il 


n?7 
(n27? + M2)? 
x {1—(-—1)"cosh M} exp { — (n?2? + M?) Dr/l?}, (6.3) 


ie 3) 
1 | (3u?-18 dud 
1 


0 


since the contribution from the ‘ Brownian velocity’ vanishes. Then 


= Poo PT. 
Rw’ | dr, | R, dr, 
~0 0 





4U4 Ke 42—1)8 =) i (—1)"cosh M}d 
= ; as my —(—1)"cosn I] 
Jo f (5 4 (n?2? + M2) ' f 
Ut ph 1 {1 3 3 | 
os 3u? — 1)8 = 4 : — : }. 
34 |, Ge PP? Me sink? a | om 
If Ul/k < 1, then _ U3la 
X3 = —__; (6.5) 
120K? 
and if Ul/K > 1, it can be shown that 
Uii3 
CF 2 8 
ae oe (6.6) 


(An estimate of X* for large values of Ul/x is given in equation (3.20) of I. The 
expressions are the same, except for the numerical coefficient.) 
It is clear that X* increases linearly with ¢ and in this sense th distribution 





becomes more skew as the time increases. However, the skewness factor X3/(X2)! 
tends to zero like t-?, and since there is no reason to doubt that the higher 
moments tend to those of a Gaussian distribution as t-? + 0 (a rigorous proof of 
this statement is difficult and will not be attempted here), we may say that the 
distribution is asymptotically Gaussian. 

A measure of the departure from normality is given by the value of the 
skewness factor. If Ul/K < O(1), the skewness factor is always small. The more 
interesting case is when Ul/k > 1. The condition that the skewness factor be 
small compared with unity is then 

Ut\t Ul F 3Ul\-3 as 

Sa sax °8 3K ania 
(This condition is similar to that given in I, equation (4.4a), but it is a little less 
restrictive.) 


7. Lateral dispersion 


The dispersion in a fixed direction perpendicular to the mean velocity can be 
calculated in a similar manner as that in the longitudinal direction. We take y- 
and z-axes perpendicular to the mean velocity, and denote by ¢ the azimuthal 
angle between the y-axis and the projection of the capillary on the yz-plane. That 
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is, if g is the velocity of fluid in a capillary, the components of velocity parallel to 
the axes are q cos@, qsin9@ cos ¢, gsin@ sing. It follows from the assumption of 
statistical isotropy that all values of ¢ (0 < ¢ < 27) are equally likely. 

The velocity of a particle in the y-direction is 


u(t) = 3U cos@ sin @ cosd + gp sin 4 cos 4, (7.1) 
and the lateral diffusivity is 
kK, = ¥ v(0) (7) dr. (7.2) 
We have, as in § 4, “ 
vv’ = 9U2 cos sin @ cos ¢ cos 6’ sin’ cos ¢’ + qpqipsin?@ cos?d. — (7.3) 


The last term on the right-hand side of (7.3) makes a contribution to x, of amount 


1 
Dsin?0 cos? 6 = -+——— | p2(1—y?) du 


3 16K 0 
_* eo 74 
3° 80 xk (7-4) 


Now to calculate 
S,(7) = 9U? cos @ sin @ cos ¢ cos’ sin’ cos ¢’, 


it is necessary to make some assumption about the statistical correlation between 
the values of ¢ in successive capillaries occupied by a marked particle. In the 
present paper, we shall assume that these are statistically independent, and also 
that the value of ¢ is independent of @. The contribution to S,(7) from particles 
which have changed capillaries in time 7 is then identically zero, and hence 


rl fl p2n 
S(7) = 90? | | H?(1— 0?) cos? d P(E 5, u, 7) d&g dud¢, (7.5) 


0/7020 





where P(&,, “,7) is given by (4.7). 
The assumption of the statistical independence of the values of ¢ in different 
capillaries is somewhat dubious when Ul/k > 1 (as discussed in more detail in I). 
Roughly speaking, this is because the lateral meanderings of the streamlines are 
restricted by continuity requirements. However, the smaller is U/l/x, the better 
will be the assumption, since the capillaries by which a marked particle leaves 
a junction are equally likely if molecular diffusion is dominant. 
(7) may be integrated with respect to 7, giving 
ton 9UF 
S,(7) = —— | lt 
8 


: ey ] oe 
2(1—y DH (57 coth M — i) du. (7.6) 


70 #0 


Then x; is the sum of (7.4) and (7.6). Altering the terms independent of U so that 


kK, = K» When U = 0 (for the same reasons as given in § 5), we have 
lau? 9022 f} 3 ee 

Ky = Kn + + | (1-7? coth M — 3) du. 7.7 

. = Oe 8 [xy al D M M?, t i 


It will be noticed that x, + x,, in accord with intuition and experience. 
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When Ul/k < 1, a UU 








1 a®?U? K s” 
kK, = 3 UL+ Si xe +0\ “3 7.9 
t 16 40 K m 3 i L ( ) 
The leading term of this expression agrees with that obtained by the random walk 
analysis. 
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, observed values (Hiby 1959). 


3 
Also shown are some measured values by Hiby (1959), in beds of spheres, for 


The theoretical values with «,, = 3« and a/l = { and ,'5 are shown on figure 2. 


large values of Ul/x. The agreement is not good but a statistical correlation 
between successive values of @, consistent with a restricted lateral displacement 
of the streamlines, could reduce the calculated values of x, considerably. The 
neglect of this correlation probably accounts for the differences between theory 
and experiment, and it is hoped to investigate this in later work. 
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Wave resistance to vertical motion in a stratified fluid 


By F. W. G. WARREN 


Imperial College, London 
(Received 14 March 1959, and in revised form 12 July 1959) 


When a body moves through a stratified fluid, i.e. one whose density decreases 
upwards, gravity waves are set up and this causes a resistance to motion. An 
axisymmetric case is considered in which a body moves steadily and vertically 
through a fluid whose density decreases exponentially upwards. The fluid is 
supposed perfect, incompressible, and unbounded in all directions. The equations 
of motion are linearized, and with a fairly general initial motion of the surrounding 
fluid, the limit of the solution as ¢t — 00 is evaluated. Transform methods are used 
to solve the equation of motion, and the methods of steepest descents and 
stationary phase are used to obtain approximate solutions. 

Streamlines and the distortion of the constant density levels for a spindle- 
shaped body are shown. The curves of resistance against a function of the 
velocity for the circular cylinder, the sphere, and a spindle-shaped body are also 
given. A criterion is given for when the maximum wave resistance for a sphere 
may be expected, and an estimate of this maximum resistance is made. 


1. Introduction 

The problem of a body rising through a stratified fluid has applications in the 
study of buoyant convection currents rising through stabie surroundings. This 
type of motion frequently occurs in the atmosphere when a thermal rises through 
a cumulus cloud. In this case the expansion of the thermal during its ascent 
causes a condensation of the water vapour within it. The release of latent heat 
warms the thermal, and this causes a further expansion and a further decrease in 
density of the rising air mass. This offsets the loss in buoyancy as the thermal 
rises to higher, less dense regions, and the upward motion is maintained (Scorer 
1958). For the purposes of this problem, the atmosphere may be regarded as an 
incompressible fluid whose density decreases upwards, and whose static stability 
parameter is the same as that of the atmosphere (Scorer 1950, 1958); and the 
thermal can be taken to have a fixed, spherical shape. What is required is a 
measure of the energy used to generate the gravity waves set up in the stable 
surroundings of the cumulus cloud as the thermal rises, so that the importance 
of this effect on the motion and growth of the thermal can be estimated. This is 
obtained in §10. 

The problem considered in this paper may also be regarded as a first approach 
to the study of the dispersion of stirring motions into gravity waves. A criterion is 
suggested for when the maximum amount of energy is used to make waves, and 
hence for when the maximum amount of energy is absorbed and dispersed by 
gravity waves. 

14 Fluid Mech. 7 
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2. The equation of motion 
Perturbation methods are applied to the basic equations of flow to obtain a 
linear equation governing the stream function of the perturbation. 
The notation employed, is, using, cylindrical co-ordinates: 
z = height (the z-axis is coincident with that of the body), 
r = horizontal co-ordinate, 
(u, w) = velocity of the fluid, 
W = steady vertical velocity of the body (positive upwards), 


: : = 1 dpy 
f = static stability parameter = — _* constant > 0, 
Po 
where /,(z) = undisturbed density, 
a* = gf, 
sa 
and k= W 
If we write 
Rs C C 
=~+u—+(w—W)—, 
Dt ct cr ( es 


Du 1ép m ; 
Dt por’ () 
Dw lép 
é 1 +--—- = —g@, 2 
une Dt * poz g (2) 
Since the fluid is incompressible 
Dp 
ae = 0 (3) 
and the equation of continuity reduces to 
C 0 
— (ru) +— (rw) = 0. (4) 
or Cz 


Equation (4) implies the existence of a stream function, (7, z,¢), such that 


low low 
uw=->~— and w=--—. 
r Cz r or 
We now write p = p,+p’ and p = p,+p’, where p’ and p’ are small perturba- 
tions of the undisturbed density and pressure, p, and po, respectively. If we 
substitute p,+p’ for p, and p,+p’ for p in equations (1), (2) and (3) and neglect 
products of uw, w, p’, p’ and their derivatives, these three equations become 


C C Op’ ; 
Pol = Wi =| nea =o (1) 
0 , oO Op , P 
Po(=,—¥ =) w+ Oz = PF (2) 
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The elimination of p’ and p’ from these equations and the substitution for u 
and w in terms of y leads to the following equation of motion: 


(;- W i (“s +2(°2)- oe (a) = 0. (5) 
ot \or\r or} dz\r oz} r Oz} r or 

The first two terms in the curly brackets represent the vorticity, the term con- 
taining / is an inertia term which arises because the density decreases upwards, 
and the last term represents the effect of gravity. 

If it assumed that a steady motion exists for which 0/¢t = 0, then equation (5) 
becomes 02 (0 (lop 0 (lew to, 12! 1) = 0 we 

sles (yo tal: zt r az|* as (Oe ait ) 

Equation (5)’ has the same form for both W and — W. Thus the equation of 
steady motion makes no distinction between ascent and descent at equal speeds. 
The assumption that there is no perturbation of the fluid at large distances from 
the body, together with the boundary condition at the surface of the body, leads 
to a solution of equation (5)’ in which a wave pattern appears below the body. 
This seems unreasonable if the body is in descent. On the other hand, if the 
assumption of no perturbation at infinity is dropped, the solution of equation (5)’ 
is indeterminate. In some cases it may be argued that the term containing # may 
be neglected because it is small compared with the other terms of equations (5) 
and (5)’. These equations then become 


[ W ) {- (-=) bo )! + a2 < (= = () (6) 
ct cz) \er\r or ez \r ez) I r or 
42 (0 (lov > sa00 1 
and : a i (-K)+=( dS Be _ +) =0 (6)’ 
022 \er \r Or cz \r oz! or\r or 


respectively. The same equations are obtained if # and g are regarded as variable 
parameters, and / — 0, g > 00, in such a way that « remains constant. 

However, the solution of the equation of steady motion (6)’ is indeterminate, 
even if the assumption of no perturbation at infinity is made. Also the cireum- 
stances in which the term in equation (5) containing # may be neglected do not 
seem obvious. For these reasons it is better to use the equation of unsteady 
motion (5), and with a general type of initial motion to look for a limiting value of 
the solution as ¢ —> 00. 


3. The boundary conditions 

If line sources and sinks are placed on the z-axis, the separating streamlines may 
be regarded, in the usual manner, as the outline of a body. The strength of these 
line sources may be written conveniently as 


—2nWf'(2)f (2) 
the accent denoting differentiation with respect to z. At points near the axis, 
Loy W ,, 
lal alae a (2)f(2) 
Fe ” 
and hence W(0,z,t) =- >the =Va say, (7) 
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where yy is defined so that the constant of integration is zero. If the axisymmetric 
body whose shape is given by r, = f(z) is a slender one, so that 7, is small, then 


wr, z,t) ” l 0 
W-—w ~ Wf(z) cz 


and W(0,2,t) if |w| <|WI. 


However, if 7, = f(z) is a boundary of the fluid, then 


UiT,...2;b) 


ee = f'(2). 


Hence equation (7) holds for a body whose shape is given approximately by 


r, =f(z), if the body is a slender one. 
The flow pattern at f = 0 is given (compatibly) as 


W(r.2,0) = Yo | 


8 
and (8) 


The initial motion described by y, and y/, is taken to be of a general nature and 
may contain waves. Alternatively, the body may be started suddenly from rest 
in a still fluid. The solution of the equation of impulsive motion (i.e. equation (5)’ 
with g = 0) is equivalent to a potential flow solution multiplied by a factor 
cos $/r. Thus the motion in the horizontal direction has a ‘wavelength’ of 2//; 
but this seems of theoretical interest only since the amplitude of the motion dies 
away rapidly in an exponential manner. 


4. The transformation of the equation of motion 
Equation (5) is transformed by means of a Fourier transform in the vertical 
direction and a Hankel-like transform in the horizontal direction: 


V(r, kt) = 


a 


v(m, k, t) dr J,(mr) ¥(r, fet). (10) 





~ 0 


The transforms corresponding to yp, (7) and y,, are 
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ee ome , Ww (2\ _ Ww Y 
Vo X%es te Xe eed TA) = F669 
respectively. If we multiply equation (5) by e“* and integrate with respect to 
z from —o to +, and assume wy vanishes at +2, we obtain 
CO ssh le flor i Seren , 0 (lc¥ 
(<+iwe) \< (-—)-2 +ipk—| +022 (-—) =0. 
ct \¢ r\r or r r . . 
If this equation is multiplied by rJ,(mr) and integrated with respect to r from 0 
to x, it becomes 


(55+ éWk) 2 ip + me) + ame] y = m(at— WE) V0). (11) 
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To derive equation (11), it has been assumed that J,(mr) (¢/ér) V(r, k, t) vanishes 
at r = Oand r = © and that Y(r, k,t) vanishes at r = 00; it then follows that 


drrJ,(mr) : | “ 
cr \ 


tos C ™ 

= = m¥'(0, k, t) —m?y(m, k, t). 
/0 PE Oe / 
Equation (11) is an ordinary differential equation with ¢ as the independent 
variable. Its solution may be written as the sum of a particular integral and a 
complementary function: y= y7,+ 7X2. Correspondingly, ‘= ‘Y,+, and 


= Wy e Wo. 
5. The particular integral 

The particular integral of equation (11) is obtained by putting ¢/ct = 0, since 
the right-hand side is independent of time. This gives y, = m‘t’,/(m” — A?), where 


k4—apks 


t= k2—k? 
| 
\ 








tan-1,/2k, 
H-7 
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FicurE 1. Path for A in the complex A-plane as k varies. 


The inverse of the transformation (10) is 
20 


Vir, k,t) =r “dm mJ,(rm) x(m, k, t). 





/0 
; ro md (rm 
Hence, V(r, k,t) = °F (6) | dm dal od 
Jo m* — A 
This integral is evaluated in Watson (1944, p. 424).+ Substituting its value, we 
obtain W'4(r,k,t) = (mi) rAH (ra) ¥,(k), 
where Im (A) > 0. 


22H) (rz) dz 
z3—)2 


gives the results immediately. 


+ The integral { taken round an infinite semicircle in the upper half-plane 
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The inverse of the transform (9) then gives 


m7 
9 /9 
2,/2 


9) 
vv, (r,z,t) = [ dk e-*? rAHP (rd) ¥ {(k), 
7 J —x 
where Im (A) > 0. The path of A, as k varies from — © to +0, is shown in figure 1. 
If £/ky < 1, this path lies approximately on the real and imaginary axes. If it is 
assumed that the body is symmetrical about its centre so that r,, = f(z) = f(—z), 
then ‘l’,(£) is an even function of k. yy, may then be written 


v= - [51m dk '¥’,(k) e***rAH(rA), (12) 


acs 70 





\ (7a if k < es 
lia, if k>k, 


where 


. kA 
and AY = gk (A, > 0), 
[A 
Ae =s5 (A, > 0) 
= k2 - v0 


6. The complementary function 
Since the particular integral is independent of time, the behaviour of yy as 
t -- «© depends solely upon the behaviour of yy, as too. The complementary 


function is obtained if we write ‘V,(/) = 0 in equation (11). This gives 


ee A,(k, m) evif + A,(k, m) eval, 


4 


where A, and A, are functions of k and m, and y, and y, are the roots of 


(y+tWhk)? (kh? —iBk +m?) + am? = 0. (13) 


Att = 0, yand ¢y/et have the given values y, and yp, respectively, and A, and A, 
are found from these conditions. We thus obtain 


A, = [Y2(XNo— X1) — Xol/[¥2— M11: 
A, = [Vi(No— X1) — Xo) /[¥1- Yo]. 
The inverse of the transforms (8) and (9) gives 


. 


, ~~ 
| dk | dm e*k mJ, (rm) LA, e%1! + A, e724]. 
x 0 


Vo = ») 

aad 
A transformation to polar co-ordinates is now made, so that one range of 

integration is finite. We write k = psin@, m = pcos@ and thus obtain 


Me ye ae 
Ws = | dé | dp p* cos 0 exp | —izpsin 6] 
\ <7 , hn 70 
x J(rp cos 9) [A,(p, 8) e%1'+ A,(p, 6) e%2'], 
ial p 
where V1.2 = —tWpsin@ +iacosé ——. 
- AN p—wsingd 
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We write 
I? se | ao | dp p* cos 0 exp [ — izp sin 6] J,(rp cos 8) A,(p, 8) e%!, 
0 


0 


, 0 (°oO 
i = — [ do ; dp p* cos # exp [ —izp sin 6] J,(rp cos @) A,(p, 8) e%'. 
Va" J —gn < 





Similarly the notations /$*) and JS are used for integrals containing A,. 

To find the limit of the integrals as ¢ > 00, we consider the p-integration first, 
and deform the paths of integration in the complex p-plane to regions where 
Re (71,2) < 0.t Two cases arise, corresponding to ascent and descent of the body. 
For ascent, it is found that if certain assumptions about the initial conditions are 
made, no poles are crossed in the deformation in the p-plane. The whole path may 
be deformed to regions where Re (y,,,) < 0. Hence in this case the integrals tend 
to zero as too. For the case of descent, the same assumptions are made, but 
poles of y, are crossed at points where y = 0, and the residues give rise to a steady 
term ast—>oo. A difficulty arises for the integrals /{-’ and JS") because it is not 
possible to deform the whole of the path in the p-plane to regions where 


Re (y;,2) < 9. 


The second integration with respect to 7 then shows that these integrals do not 
have limits as t-> 0. However, if the assumption is made that //k) < 1 this 
difficulty is partly overcome in the sense that the integrals /(-) and J") taken 
along the deformed paths remain negligible for a certain interval of time. During 
this time the integrals are approximately equal to the residue terms from the poles 
crossed by the deformations. From these residues the limit of yy, is obtained. Thus 
it may be shown that 

| 0, if W>0 | 
ok 


lim y= (14) 


ran |J2m [de WoC) PA (rA) sin ke if W<0. 


The remainder of this section gives a proof of the expression (14) and discusses 
some points which arise from the preceding paragraph. 


Case I. Ascent W > 0 

It is assumed that A, and A, are O(exp[b|Imp|]p-*-4) as |p| > 00, where 
A> 0Oandb >0. Equation (13) shows that if y = 0, then m*? = A. Hence y, has 
poles at points where y = 0, A, has a pole at y, = 0 and A, has a pole at y, = 0. 
Figure 2 shows the p-plane. Regions where Re (y) < 0 are shaded, and the poles 
of A, . at y;.. = 0 are shown at the points marked I’. y has a pole at p = ifsin@ 
and a branch point at the origin. A cut is made along the imaginary p-axis to 
join these points, so that y is single valued in the cut plane. 

We consider the upper bounds of the integrands of the integrals /{*, ete., taken 
along the deformed paths. Upper bounds for the integrals taken along the 
deformed paths may then be found. The range of the 4-integration is divided into 
the ranges (0, +¢) and (+¢, + 47), where € is a small positive number. 

+ The transformation p — p(y) is not a simple one and it is better to keep to the inde- 
pendent variables (p, @). 
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Discussion of the integral I? 
(i) No deformation is made in the interval 0 < 4 < €, and it may be shown that 
Poe 
dp < M, for some fixed M. 


70 


(ii) Iné < @ < $m, the path of integration in the p-plane is deformed as shown 


in figure 3. 
18 sin@ VR 
Cut ry, =0 


C/ pi B sind 
I" " I; 
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FiGuRE 3. Case I. Deformation for Jj ine < 4 < 4n. 


It follows that 


dp < Ke 
/ R,Rex 
by suitable choice of #,, and that 
! a K 
dp < = : 
JOR; Wtsin 0 — 2(r + |z| +b) 


if Wtsin@ > 2(r+|z| +5). where K is a fixed number. 
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Discussion of the integral I~’. 


(i) In the interval —eé < # < 0 the path is deformed as shown in figure 4. 


[t may be shown that dp, < M eo'r+='+, where c is some constant which 
20,0 
depends on f# and k, only. 
(ii) In —47 <6 < —e, the deformation is similar to that used in /{*, see 


figure 5. 





Via 








FicureE 4. Case I. Deformation for I{’ in —e < 0 < 0. 


R, 











FicurE 5. Case I. Deformation for I? in —47 < 6 < —e. 


In this case it may be shown that 
: ’ K 
dp| < Ke+ __. ‘ 
J ori Rox Wt sin 0 — 2(r + |z| +6) 

Similar results hold for the integrals /$*’ and JS’, when the paths are deformed 
in a similar manner. If we collect these results for the integrals /{*’, ete., and 
consider upper bounds for varying 7, we obtain the following upper limit for the 
d-integration: 


ehn * 


nwo —_ 


dé dp < Ki +exp [ck (r+ |z| +5)]) 
J —tn ~ paths 
where Wt,sine > 2(r+ |z| +5). 

For poles crossed by the deformations there will be contributions from the 


> 


: | 
let Wt, sin ¢€— 2(r + |z| +6)} 








residues, which give terms like 


be 
[ dO © evrt. (15) 
/-4n 
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Here © is some function of @ and a pole occurs at p = p,(9), say, where 


Y = Y(Pp9) = ¥p(9). 
If Re (y,) < 0, the expression (15) will decrease with time, but if Re(y,) > 0 (e.g. 
for the deformation used for /5~), this may not be so. In general, if it is possible to 
deform the path of integration in the @-plane to regions where Re (y,) < 0 (where 
Y, is such that Re (y,( + 47)) < 0), the expression (15) will decrease with time. 
However, this deformation is not always possible, and, for example, we may 
have Im(y,) = constant = WI, say, and Re(y,) = 0. There is then a steady 
fluctuation present throughout the motion. This type of motion is easier to see if 
we use equation (6). The assumption that 
b oc et xJ,(nr) eo" 
then gives (L+k)? (k? + m?) = k2m?. (13’) 


Hence equation (6) has a solution 
(PaO 
=reiWt | dk (k) e-** mJ,(rm), 
70 


— 


where m satisfies equation (13’) and J is fixed. 

This fluctuating type of motion does not seem to be of practical interest since it 
must be present in the initial conditions and it is not easy to see how this motion 
could arise in the first instance, unless a fluctuating type of disturbance were 
present ° 

With regard to the /-integration for the residues in the case where 

Im (y,) = constant, Re(y,) > 9, 


the fluctuation would increase with time. In this case the linear theory is inade- 
quate to say what eventually happens. 

If it is now assumed that poles of the above description do not occur in A, and 
A, and hence in yy) and Yo, i.e. that certain types of waves are absent from the 
initial condition, then 


= 


= “olf cciedh ] | 
+ eX KC P+ 12 ij Re 
K (1 exp [cko(r + |z| +)])¢ aay ear ew 
| 


This holds for all (7,,z,,¢) if 7, +|2z,| < r+ |z| andt >t). Hence y,> 0 ast—>ooif 


W > 0 for all regions of the flow. This result holds for any # > 0. 


— 


Case II. Descent W < 0 


The same methods are used as in Case I, but it is no longer possible to deform the 
paths of integrations for the integrals /\~ and /$") in the p-plane so that they lie 
completely in regions where Re(y) < 0. The singularities of y prevent this. 
Figure 6, which corresponds to figure 2 of Case I, shows the p-plane. 

The paths for the integrals /\-’ and J5") are deformed as shown in figures 8, 9 
and 10. In the regions were Re(y) > 0, the paths are deformed so as to pass 
through a col as in the method of steepest descent. An upper bound for the 
contributions from segments of the paths where Re (y) > 0 is then made. The 
track of the cols for varying @ is shown in figure 7. 
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The position of the cols is given by p = p,(@), say, and at these points 


p V(Pes ) = y(9). 


: j 3,/3 B 
We write cosec #* cot 6* = — o (47 > 6* > 0) 
’ 0 
; 1p 
and cosec 0** cot 0** = 5 i. (47 > O** > 0). 


one 
Then if |6| > 0*, it may be shown that two cols lie on the cut; and that if |@| > 0** 
both zeros of y lie on the cut. We suppose that / and g are variable parameters 

















FicurE 6. Case II. Paths of Re(y) = 0 in the complex p-plane. 








6=0 [+ (1Bky) }: 





FicurE 7. Case II. Path of the cols for varying @ in the complex p-plane. 


with x, k, and W held fixed, and that 0 < £/k, < 1. The range of the 0-integration 
is divided into the ranges (0, +¢€); (+¢€+4a746); and (+4746, + 4m); where ¢€ 
and 6 are small positive numbers. / is chosen so that 47—6* < é in order to 
restrict the contributions from the paths in regions where Re (y) > 0. 


Discussion of the integral IS*? 


(i) There is no deformation of the path in 0 < @ < € and it may be shown that 


hy 


| dp < M. 


~0 
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(ii) Ine < @ < 3m, the path is deformed to Re (y) = 0. On this path it may be 
shown that dy|dp 
“sa > 
W sind, 
With the assumption that A(p, 4) is of bounded variation on any finite length of 
the path, it follows that 
i K exp [cky(r + |z| +5)] 
dp < = 
Wtsiné 


fa C 
ff KID LL 


x C 
Im (y) = Im (y,) 


a: Me. 


Similar results hold for J[~’. 











FicurE 8. Case II. Deformation for I$ in 0 < 6 < e. 


Discussion of the integral IS*? 
(a) Contributions from the paths which lie in regions where Re (y) < 0 
(i) In 0 < @ < €, the path is deformed as shown in figure 8. 
It may be shown that 


| dp| < Mexp [cky(r + |z| +5)]. 
/OX+Yo 


(ii) Ine < @ < 42-6, the path is deformed as in figure 9. Here 
r | Kexp|[ck,(r+ |z|+6 
| dp| < J | “ol? al ) 
Jox | Wtsin@ 
as may be shown by a discussion similar to that used above for the integral 
I*?, (ii). 
If £/ky < sind, it may be shown that 
dy/dp ] 
*S > 
Wsinf|” 2 
and that the integrand is uniformly bounded with respect to # and @ on the 
deformed path. 


It then follows that 
if | 2 K exp [ck,(r + |z| +5)] 


d —_ 
Jrz v Wtsin@ 


where |OZ| = kg, say. 
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As in Case I, integral /{*), it may also be shown that 


; | , K 
dp a 
Jaa! | ~ Wtsind—2(r+ |z| +5) 


ip 
| 


ae Me. 


(ili) In 37-6 < 4 < 4m, the paths are deformed in the three subintervals 


ae ec OaO Oc Oc Om OP <<; 


> 


/ 


I-39) 


Vee i 
Cf: Im (y) = Im (y, 
x, 


i} 











+) 


Figure 9. Case Il. Deformation for IS" ine < 0 < 37. 
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:1-3<9<8 §* <9<9** O° <9<in 
Figure 10. Case II. Various deformations for I; in 31-6 < 6 < 3a. 


as shown in figure 10. For points on the segment Yoo of the deformed path. 
where |OY| = ky, it may be shown that 
dy/dp Pe J 


—- a 
Wsinf| 2 


Similar arguments hold for | dp| as in (ii) above. 
JV | 
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(b) Contributions from the paths which lie in regions where Re (y) > 0 

If we consider the path of the col, we may show that in the interval 

0<0<6*, Re(y,/W) < 5\(fky); 

and that in 0*<A<4n, Re(y,/W) < $2. 

(i) In 0 < 4 < 4n—6, it may be shown that the length of the path XY is 
less than K J (Bk). 

(ii) In }1—6 < 6 < 4m, the length of the path in the region where 

Re(y) > 0 is O(kp). 

It may also be shown that the integrand is uniformly bounded with respect to / 
and 4 on the deformed paths. From these results it may then be shown that 


ohn . 
| dé | dp < K{V(Pko) +6] exp [cky(r + |z| +b)] exp [5W ./(Pk,) t]. 
0 7 XY 


The integral /{-) may be discussed in a manner similar to I$*?. 
From the above results for the integrals /{*, ete., we obtain the following upper 
limit for the 7-integration: 
. | Ip| < Kexp[cko(r + |2| +5)]| : 
( dp| < Kexp(ck,(r+ |z| +6)]\e+ — —-+ 2. 
f PLNo ra {| Wtsine Wtsine—2(r+ |z| +5) 


J —4n / paths 


+ [V(Bko) + 8] exp [5Wy(Bky) t } ; 
This holds for all 0 < £ < #, for some Pf > 0. (c does not increase if # decreases.) 
[t also holds for all (7,,2,,¢), where r, + |z,| < r+ |z|, and t > fo. 

If we choose t, >t) and # such that J(fk)) < L/S5Wt, (i.e. V(A/ky) < 1/52t,), 
then iis ‘ 
~ dé dp| < Kexp|[ck,(r + |z| + 5)] (3¢ + 46) (16) 


J—4n / paths 





during some interval of time (t,, t,), f: > t. If (e, 6) is chosen so that the expression 
on the right-hand side of (16) is small compared with the average value of 


|v, + residue terms|t 
in the region r, + |z,| < r+|z]|, then 


' rh 
| dé | 


jv —$7 ~ paths 


. 


dp| 


may be neglected in the interval t, < t < ft. This interval may be called a ‘ quasi- 
7 (20 
dé | dp for large t may be found from 


/ —$7 /0 





limit’ interval. An estimate of yy, = 





the asymptotic expression of | dp. The method of steepest descent through 
x¥ 


the col gives the first n-terms in the formt 


n [hn 
y 
ae 


| dO O, exp [y,(9) tt-7-8, 
r=1J—}n 


+ yw, and the residue terms change little as # > 0. 
t A modification is needed at 6 = 6* since dy/dp = 0 there. 
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where the 0, are functions of #, and Re (y,( + 47)) = Re (y, (0)) = 0. A different 
form of deformation of the path in the p-plane occurs at 0 = 6* and hence not all 
the Q, are holomorphic at this point. Hence complete deformation of the path of 
integration in the #-plane to regions where Re (y,) < 0 is not possible,+ and 





gives terms 
exp[}(,| =) 


o( me 


where r is some integer. Hence |1/,| diverges as > «. and eventually the motion 


increases without limit. 








! 

up + Ati 

Pe = atl at Baie 
g=0 7] 6 one a 
—S= J=0 
whe 6<0 8>0 ho <0” ns 
— 
a (b) 


FIGURE 11. Paths of integration for k, (a), and A, (6). 


The assumption that second-order terms are small no longer holds at this stage 
and the linear theory seems inadequate. Physically, the medium is a dispersive 
one in which local disturbances are propagated outwards in all directions. The 
gradually increasing inertia forces retard the propagation to deeper levels, and 
waves ahead cannot disperse so rapidly in descent as in ascent. 

Similar assumptions concerning the poles of vy, and ¥, are made as in Case I. In 
the case where the motion is started suddenly from rest, a pole of yj occurs at 
p = ifsin#. It seems better to treat this pole separately, since the arguments used 
earlier would not hold in this case. 

With the above assumptions, the limit of the complementary function in the 
‘quasi-limit’ interval is found from the residues of the poles of x, at y = 0, which 
are crossed by the deformations used for /{-’ and J,"). If we substitute for the 
residues, after a little manipulation we obtain 


-. ie j ‘kez 
- IE ir | dk'¥’(k) e-* AJ, (rA) sgn 4, 
« JL 
where L is the path in the k-plane shown in figure 11. Figure 11 also shows the 
corresponding path for A. 
If fr < land fz < 1, then these paths may be taken to lie approximately on the 
real axes. With the assumption that ‘’,(k) is even, this leads to the expression 


Ko 
J2n [ dk V(k) rAJ,(rA) sin kz. (17) 
Jo 


+ Other paths may be chosen, but a different form of deformation occurs at some stage 
because the contours of y change their form at 0 = 9* and at 6 = 6**. 
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In conclusion, for Case II, if £/k, < 1, there is a region around the body where, for 
a given time interval, yf, is given by the expression (17). The length of this interval 
becomes infinite as //k,y—> 0. 


7. The complete solution 
The results of §§5 and 6 give the following expression for y as t > : 


yf = : Im | dk F(k) rAHP(rA) exp [7 (sgn W) kz], (18) 


i 
sa 70 
a 


where Fi(k) = | dz cos kz { f(z)}. 
/0 

Hence the wave pattern produced by a descending body is the same as that 
produced by an ascending one. The waves are downstream in each case. When the 
initial conditions are those of a steady, numerically greater velocity, and yp is 
given by expression (18), it may be shown that the same wave pattern eventually 
emerges. 

Methods similar to those used in §§ 4-6, may be employed to solve equation (6). 
There is no essential singularity of the integrand in the corresponding double 
integral expression for i/,. The paths of Re (y) = 0 coincide with the real p-axis. 
The poles of y, also occur on the real p-axis, and the double integral expressions 
for yy, and yy, are improper integrals. If principal values of the integrals are taken 
we obtain after one integration 

; i 
dk F(k) rAH{P(rA) cos kz 


#0 


iW 
f, = Im 
») 





‘ ie 
and lim yy, = . Im | ‘dk F(k) rAHP(rA) i(sgn W) sin kz. 

t>a - 70 
These two results differ from the corresponding expressions when / + 0. However, 
the sum yw = y,+7/, is the same in both cases, and in this sense the limit of the 
solution of equation (5) as / > 0 is the same as the solution of the limiting form of 
equation (5) as / —> 0. 


8. Asymptotic approximation to the motion 

An asymptotic expression for the streamlines and the distortion of the constant 
density levels can be obtained by the method of steepest descents. The displace- 
ment of the streamlines is given approximately by 


- ' 
u 
4 i 
A,= | dz—, = —-~, 
ae F VW Wr 
and the displacement A, of the constant density levels by 


A, = | dex. 


“2 


If we write k = k,sin € and use the asymptotic expression for H{)(2) 


H'(x) = [pcexplite- i77)). 
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then expression (18) yields 
k . 
a, & - ff im | dé F(k,sin &) sin £ cost £ exp [tky sin £(z +r tan £) — 377i], 
7 ~M 
when W > 0. M is the path shown in figure 12. It is deformed to pass through the 
col at C. This col contributes an exponentially small term to the asymptotic 
approximation to A, and is omitted. It may be regarded as the potential flow 


displacement. 








Saddle point 
vo 





FicurE 12. Path of integration for &. 


When z < 0 a saddle point appears on the real axis between the origin and 
(477, 0). The method of stationary phase then gives 
k 


“0 x» —- () 


sin £* cos? £* cos 


A, = Lk F(kysin £*) 


hk ate 
zk, sin & ) if 
* 7 a{(2+8in? &*) 


‘ 2 0% 
1+ cos? é 


Oo fF z>0 
for the asymptotic value of A.. Here &* satisfies the equation dG/dg = 0, where 


G(€) = irkysin g(tang+-), 
and where z/r = tan ¢, say. This gives 

tan? £* + 2tan é*+2/r = 0 
or sin 2£* + 2 tan (€* —¢) = 0. 
A similar result holds for A, : 


4 1k, F(k,sin &*) 
po» (2+ sin? £*) 
v= , > J 


sin? £* cos £* cos (rk, tan? €* sec E*) if 2 < 0 


xa 1 2> 6 


If W < 0 the reverse holds with A, and A,, asymptotically zero if z < 0. 
15 Fluid Mech. 7 
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FIGURE 13. 








Svreamlines for a spindle-shaped body. 
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If €* is used as a parameter we may write the displacements as 
1 ; , , ' 
— x amplitude factor x cos {(r or z) x phase factor}. 
; ; 


The results for a spindle-shaped body whose shape is given by 
x 

acos— if |z| <b 
2b 

=0 if |z|>b 


II 


J () 


have been plotted and are shown in the figures 13 and 14. a and b are the half- 
breadth and half-length of the body and kj* has been used as a unit length to 
measure distances from the body. ak, has been used as a unit displacement to 
draw the streamlines and the constant density levels. It may be shown that the 
slopes of the curves joining the crests and troughs (shown dotted in the figures) 
are asymptotically — cot £* so that the approximate rate at which the disturbance 
spreads outwards is }rW. A value of $7 for bk, was used in the calculations. 


9. The wave resistance 

The wave resistance to the motion of the body is found from a consideration of 
the pressures on its surface. If the static thrust is neglected the upthrust is 

" rb ‘ ds : 
R=2n| dzpsinn— f(z), 
‘= dz 

as shown in figure 15. p is the perturbed pressure. Euler’s equation of motion 
relative to the body for horizontal motion under steady conditions is 


when products of small quantities have been dropped. The assumption that 
1/p = 1/p> gives 
, [2dret Wp (° 9 rem an iin 
p=We| “= | dk k2F(k) H2(rA) exp [i(sgn W) kz]. 
Jr 1 Ce =“ Jo 

Further, we assume that the body is slender so that ds/dz ~ 1 and siny ~ f’(z). 
If we write Po(z) ~ po(1—/z), where py = density at z = 0, then we obtain 

R lg , er ‘eo oP . (1) 
= Im | dz(1 — fz) f(z) f' (2) dk k?F(k) exp |i(sgn W) kz] Hg? (r, A). 


7) ~0 





7),W? 
Since f(z) is even, this becomes 


R 
27p,W2 


*Ko rh 
= (sgn W) | dk k*?F (kk) | dz f'(z) f(z) sin kz Jo(7,,A,) 


0 


j 
/9 ny 


D Cb 
~pIm | dk k?F(ke) | dzf'(2) f(z) cos kz HP(rq A). 


/0 /0 


a 


The latter integral represents the effect of the change in virtual mass, and is 
15- 


9 
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neglected because it is multiplied by /. If 7,A, is small over a large part of the 
range 0 < k < ky, s that J(r,A,) = 1, then approximately 
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FicureE 16. Curves of R, against v. 
Wave resistance 


fa'g po 





_* 
R,= 





FicurE 15. Calculation of bie bV(9P) 


the wave resistance. mW 


If the shape of the body is written as r,, = af,(z/b), where a = half-breadth and 
b = half-width of the body, the wave resistance is given by 

R ; FY oie 

p= aa | dk RPh}? = Ry, say, (19) 
Po® OP mV" Jo 


where v = bk,/7 and 


Fy(k) = | dzcos ket f,(2)}*. 
#0 
For small yp, R, oc v®. In many cases the right-hand side of equation (19) will tend 
to zero as y-> #, and hence, in general there will be some condition of maximum 
resistance. 

Wave resistance curves for the spindle-shaped body, the sphere and the 
circular cylinder (whose shape is given by f(z) = a, if |z| < b, f(z) = 0 if |z| > b), 
are shown plotted against v, that is against a function of velocity, in figure 16. 
Many assumptions break down for the circular cylinder and its resistance curve 


is of theoretical interest only. 


10. Applications and conclusions 


The rise of a thermal through the atmosphere has a velocity given, approxi- 


mately, by Scorer (1957) W = (Ba), 
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where B = buoyancy = density difference/density, and 2a = width of the 
thermal. To estimate the wave resistance it is assumed that the thermal is a 
sphere whose diameter is the width of the thermal. The ratio of the wave resistance 
to the buoyancy is then 


resistance 3 af ap _(ap\? 
= Ral — = 0-05( D> 3 
buoyancy 476 AN 7B, B, 
roughly, ifaf/B < 3. For the following values of the constants (which may occur 


in the atmosphere) 
ao=1l@cm., B=10*, £= 10% cor, 


this ratio is about 0:0005. Thus in these cases the wave resistance to motion is 
negligible. This remains true over a large range of values of B and #. On the 
other hand, the following values of the constants: 


a=10cm, B=3.10%, ~= 10’ cm, 


give a ratio of about 0-5. 
As a rough guide, a sphere will experience most wave resistance at a speed 


given by : 
; av gp ms 
_—- = . 


Vv 


The ratio of the wave resistance to the weight of the displaced fluid is then 
about jaf, where a is the radius of the sphere. 

Hence in general the dispersion of energy due to wave motion is small. How- 
ever, when the situation is specially favourable to the surrounding fluid, and the 
forces involved are small, this dispersion may play a major part in the behaviour 
of the fluid, as in the case of the larger thermal with a smaller buoyancy (e.g. 
towards the end of the ascent). 

In a similar manner, it may also be possible for these waves to affect the 
behaviour of turbulent eddies. For instance, if the eddies are large, and if their 
outer regions move slowly, a considerable proportion of energy may be dispersed 
by gravity waves. A criterion for this energy disperion to be a maximum would be 


ay gp 
_ 


where the number » depends on the shape of the eddy, a is a typical dimension, 
say the eddy size, and W a velocity typical of its outer regions. 


I would like to express my thanks to Dr R. 8. Scorer for the very many helpful 
discussions and guidance he has given me throughout this work. 
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The viscosity of a suspension of spheres 


By A. D. MAUDE 


City of Liverpool College of Technology 
(Received 6 May 1959) 


[t is shown that in Stokes’s flow the perturbation field, due to the addition of one 
more sphere to a shear flow of a fluid containing a number of non-interacting 
spheres, has the property that the total additional shearing force, acting on any 
plane normal to the direction of velocity change, is zero. However, the perturba- 
tion velocity, integrated over such a plane, takes a constant value, positive if the 
plane lies on one side of the sphere and negative if it lies on the other side. It 
follows that the effect of all the spheres is not to alter the shearing stress at all, 
but to reduce the mean shear by a factor 1 — 2-5c, where c is the concentration. 
This suggests that Einstein’s viscosity law should be altered to 7 = 7)/(1 — 2-5c) 


when c is not small. 


1. Introduction 

If a number of spheres is added to a fluid the viscosity of the fluid appears to 
increase. Einstein (1906, 1911) has shown that if the density of the fluid is the 
same as that of the spheres, if the Reynolds number is sufficiently low for all 
inertial effects to be neglected, if the spheres are much smaller than the vessel 
containing the liquid, and if the spheres do not exert an attractive or a repulsive 
force between themselves, then for very low concentrations 


N= No(1 +250), (1) 


where 7 and 7, are the viscosities of the suspension and fluid respectively, and 
cis the volume of spheres in unit volume of suspension. Many attempts have been 
made to extend this equation to higher powers of c for the case of equi-sized 
spheres. In most previous work, simplifying assumptions have had to be made, 
and the present work is an attempt to extend (1) without the use of such 


assumptions. 


2. General considerations 

The liquid is considered as being sheared in an idealized viscometer consisting 
of two infinite planes y = 0 and y = H, the latter being moved with a velocity 
U, in the positive direction of x by a force F per unit area. The viscosity is then 


ny = FHIU,. 


given by 


The introduction of spheres into this liquid will disturb the flow pattern, and it is 
convenient to consider the disturbance such that it vanishes at infinity. It is 
shown below that this implies that the average value of F is unaffected, although 
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the average value of U, is reduced. The perturbation velocity in the x-direction, 
v,, will cause the velocity to vary over any plane y = const., but when the plane 
considered is far from the spheres, and when the perturbations of many spheres 
are superimposed, the variation over the plane will vanish. This variation over 
the planes y = 0 and y = H is called a wall effect, and will be neglected in the 
present work. 

The flow pattern may be calculated in theory by the method of successive 
reflexions. The first-order perturbation velocity caused by each sphere placed in 
the original flow is first calculated. The sum of all these perturbation velocities is 
then taken as the initial velocity field, and the second perturbation field due to 
each sphere placed in it is calculated. Third and higher-order perturbation 
velocity fields can be calculated similarly, leading to an infinite series of terms. 
It will be shown below that second and higher-order perturbation velocity fields 
cannot affect the viscosity of the suspension. 

It appears to be difficult to prove the convergence of the above series, and, 
while it is expected to hold for low concentrations, there may well be a critical 
concentration above which it diverges. 

Lamb (1945) gives a general solution of the linearized Navier Stokes equation 
as follows: 

; ij #* dG, apt 0 p, \ 
ey \2(2n4+1) Ca (n+ 1)(2n4 1) (204 8) Oa r2"+1| 


n\ 9 
| (2) 


with two similar expressions for v, and v, obtained by cyclic permutation of 
x,y,z. Here r? = 2?+y?+2?, and p,, ¢,, X, are arbitrary solid harmonics of 
degree n, and the pressure p = Xp. 

If the origin is now taken at the centre of the added sphere, the restriction that 
the perturbation velocity and pressure must vanish at infinity implies that n < 0. 
For mathematical rigour, a proof is needed showing that it is always possible to 
satisfy the required boundary condition at the sphere surface; such a proof could 
easily be constructed using the method of §7, but will not be given here. 


3. Lemma 1 

The force in the x-direction on any plane y = const. + 0 due to a perturbation 
velocity field is zero. 

The force on a plane y = const. will be given by 


Poo Proc Cv Ov 
y) | | (< M4 ) ded; 
J 2= oe EE L 


oa 
ox oy 


and as this must be independent of y for there to be no accelerating force on an 
infinite slab between two planes y = const., only harmonics which give rise to 


v, and v, of degree — 1 need be considered. 

Neither ¢ nor X can give rise to terms in v of degree — 1, and if the sphere moves 
so that it experiences no resultant force, the significant term in p,—that is, p_.— 
must be zero. 
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4. Lemma 2 

The value of | |v,dxdz in a perturbation field is constant for all positive and 
negative values of y. 

From lemma | the total force in the a-direction on any plane y = const. is zero. 


That is, oe om P 7 
Cov ov 
| | (G4+ 52) axde = 0. 
Jz=-wJr=-w \OU OY, 
Pa Px ov rx 
Sut y) | —“drdz = 9 | [v,]e--— 0. d2 = 0. 
wt 0 of Zt 2 Ol JZ L 
Po lie ov 0 ex 
Therefore | | —"dadz = = | | v,.dxdz = 0, 
Jt OHO wJT  C Y Cy JZ OJ x 
f* ax > ox 
and | | v,dxdz = const. (3) 


5. Lemma 3 


The average value of (cv, /0x + ¢v,/ey) and of 02p/ca cy over the surface of any 
sphere r = const. is zero. 

If expressions for 0v,/¢.x, 6v,/Cy and ¢*p/cx cy are found from equations (2), and 
are integrated over the surface of a sphere, only a very limited number of values 
of n will be found to give a non-zero value, and in each case n > 0, so that this 
term will not occur in the perturbation velocity. 

For example, consider the contribution of ¢,, to the integral of ¢v,/cy over the 
surface of a sphere, which is ae 
( “On 
— dA. 

CUCY 

As @,, isasolid harmonic, ¢7¢, /¢2 Cy is also a solid harmonic, and by the orthogonal 
properties of harmonic functions this integral will be zero unless 07¢,,/Cx cy is of 
degree 0 or — 1. Therefore, v,, = ¢¢,,/¢x will be of degree 1 (the degree zero not 
giving 07, /cxcy of degree — 1). Butasv, must vanish at infinity this term cannot 
be included in the perturbation velocity field, so that no velocity due to the 
term ¢, can contribute to the integral. 

A similar argument applied to each possible harmonic proves the proposition. 


6. Expression of the significant parts of the integral (3) in terms of 
Pn Pn» X, 
As the integral (3) is a constant, only terms in v, of degree — 2 can contribute 


to it, so that we are restricted to p_;, 6_,, X_,. Let 
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If we express p_;, 6_,, X_, in spherical polar co-ordinates, using the axis of y as 
the line 0 = 0 and the positive (y, z)-plane as the plane w = 0, considerations of 
symmetry show that the only terms which contribute to the integral and give 
a value which is different for positive and negative values of y are as follows. 
rs ae ee p_3 = Arcos wP3}(cos 4), 


where P3(cos #) denotes the associated Legendre polynomial cos @ sin @. This gives 





(ox 0 (7A/3y for positive y, 
| v,,dxdz = ; . ; 
Je=—0 Jz=-a |—7A 3y for negative y. 
Secondly, there is X_, = Br sinw P}(cosw), 
which gives 
ox pax fn 27B for positive y, 
| | v,,dxdz = ; : 
Je=—0 Jr=-« \27B for negative y. 


[t should be noted that terms which make the integral give the same value for 
positive and negative y may be neglected since they do not affect the rate of shear, 
and hence the viscosity, of the suspension. 


7. Expression of A and # in terms of the original flow 
To avoid the complication of the movement of the sphere we may impress a 
translational and a rotational velocity upon our axes, and thus upon the original 


flow, such that 
[47a curl v],_, = 90, 


| 67yav + 7a grad p],_» = 0. 


The bracket [ ] is used to indicate that the function is to be calculated for the 
original flow, and the suffix 7 = 0 to indicate that its value at the origin is to be 
taken. Here a is the radius of the sphere. 

By Faxen’s theorem (proved simply by Pérés in 1929), this implies the absence 
of movement or rotation of the sphere if no external force or couple is exerted 
on it. This choice of axes will not affect the perturbation velocity field. 

Lamb gives e 

nr? 
. by Vand 
TU, = — 2 — p, + UNQ,, (4) 
7 9 2(2n +3)! ” sii 
re ea: \ ‘4 5 
rceurl,v = Xn(n+1)X,. 5) 


The left-hand side of both these is determined by the original flow at r = a. 
In both the original flow and in the perturbation flow 
i l © (y2 1 Fe in 0) 1 dv, . 
aiv ¥ = —— (f°v,) + —— — (0,8 + Ey ys oS 
ror’ ™  rsin@cd: ° rsin@ cw 
As v, and v, are equal and opposite in the original flow and in the perturbation 
flow at r = a, 0(r?v,)/Or must also be equal and opposite in the two flows at r = a. 
But from (4) we have 
0 1 n(n+3)r? 


a (rv) = * 2(2n +3) Pyntun(n+1)¢,, (6) 


and the left-hand side of this is determined by the original flow at r = a. 
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If both sides of (4) and (5) are multiplied by cos wP}(cos @) and are integrated 
over the surface of a sphere, the orthogonal property of surface harmonics allows 
us to write 








Pa 2a 
L=- | | [rv,],-_. CoS WP3(cos 7) dw d cos 6 
J-1/0 
671A 3670 
~ Ba 5 a?’ 
"1 ena 
I, =— | — rv, cos w P}(cos 4) dwd cos 4 
7-1/0 or r=a 
72nC 
~ § @®’ 


where C is the constant of the harmonic 


¢_3 = Cr cos wP3(cos #). 
Similarly, 
re 21 
i -| | [rcurl, v],_,, sinwP}(cos@)dwd cos @ 
é ba? 
87 B 
3 ae 
Fi aco 22 5 na 303 
a Chis gives A= ! (24,+1,), B= ies 
12 7 sj Sim * 


As I,, I, and J, are determined by the original flow, A and B may be found for 
an arbitrary original flow. A similar procedure could be made to yield the constant 
for each of the harmonies in the perturbation velocity field, and thus to solve the 
problem of a sphere in an arbitrary velocity field. 


8. Expression of /, and /, in terms of the original flow at the origin 
Let the original flow be expressed as a series of spherical harmonics in the 
manner of equation (2). As it contains no singularities within r = a, only 
harmonics with n > 0 are included. 
Consider the harmonies 
[po] = [A’r? cos wP3(cos 4)], 
[X,] = [B’rsinwP}(cos4)], 


! 


[J] = [C’r? cos wP}(cos A)]. 


Then in a similar manner to the last section we get 


127 a4 247 , 
L==__—s’-—av’, 
35 7 5 
12ra* ., 727 ; 
L,=--— A’ — a?C 
a ‘YY oD 
S7 
ol 
| # = 3 ab’, 
Therefore A = —a°A’— 10a37C", 


and B= —-aB’. 
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Consider now [¢7/p/0x Cy],,_ ». Clearly, only harmonics of degree 2 can contribute 
to its value, as higher degrees will vanish when r is put as zero. Application of the 
double differentiation to all possible harmonics of degree 2 shows that 


ro Ov, ov 

similarly, Ba BE = 2C", 
. COX CY Ir=0 

and [curl vL_g = B’, 


which is zero by the hypothesis at the beginning of §7. 
It follows that the difference between the values of | |v,,dadz over each visco- 
meter plate is altered by an amount 


9. Viscosity of a suspension 

If there are n spheres in unit volume, the value of U,, is reduced by an amount 
nHS, where S is the average of S averaging over all the sphere centres. But from 
lemma 3, the perturbation velocity cannot affect the average value of S, so that 
it is determined solely by the original flow. That is, 


Therefore U.. is reduced to 
Pr 


U (1 —3n $na3) = U,(1 — 2-5c). 
The viscosity is thus 
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10. Limitations of the theory 

[t has been implicity assumed in the last section that no ordering has taken 
place. That is, any one sphere will be situated in perturbation fields which arise 
evenly throughout space. In fact, some ordering will occur. It is well known that 
the probability of finding a second sphere centre a distance r from a given sphere 
centre is zero up to r = 2a, rises to a maximum and falls again, giving an oscil- 
latory function of r of rapidly decreasing amplitude. The maxima become sharper 
as the concentration increases. However, except near to a wall, no ordering of 
direction occurs, and as the average values of ov,/éx+0v,/¢y and 0?p/éx cy are 
zero for any value of r, the ordering with respect to distance will not affect the 
result. 

The possible divergence of the series mentioned in §2 might produce a change 
of the law at some concentration. Experiments by Higginbotham, Oliver & Ward 
(1958) show that the viscosity of a suspension of spheres obeys the law 


1 = No/(1— Ke) 
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up to a concentration of 28°, above which concentration there is a marked 
departure from this relation. The value of K obtained was slightly less than 2-5, 
It is suggested that this change of the law at c = 0-28 is due to the divergence of 
the series at higher concentrations. 

Equation (7) above agrees with Einstein’s result (equation (1)) if the square 
and higher powers of c are neglected. It also agrees with a theoretical expression 
derived by Kynch (1956) on what he calls a ‘rigid envelope’ model. Kynch’s 


expression is m ee 
N = No( 1 + 2:5c + 6-25c2 + ...). 

However, equation (7) does not agree with most of the other theoretical expres- 

sions which have been published. 


The author wishes to express his thanks to his wife and to Mr H. Mulholland for 
help in preparation of the manuscript of this paper. 
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Force-momentum fields in a dual-jet flow 


By DAVID R. MILLER* AND EDWARD W. COMINGS?+ 


Purdue University, Lafayette, Indiana 
(Received 17 April 1959) 


Measurements of mean velocity, mean flow direction, normal turbulent stress 
in the direction of flow, and mean static pressure are reported for the subsonic 
flow field generated by identical twin jets of air issuing from parallel slot nozzles 
in a common wall and mixing turbulently with ambient room air. At the low 
nozzle velocity employed (72 ft./sec), the two-dimensional plano-symmetric flow 
was effectively incompressible. Since the end walls prevented interjet air entrain- 
ment from the surroundings, a region of highly convergent flow was formed near 
the nozzles. In this region, contour maps clearly reveal (1) the sub-atmospheric 
static pressure trough that accounts for the jet convergence, (2) a free stagnation 
point on the plane of symmetry, (3) stable symmetrical contrarotary vortices 
which recycle air on the concave side of each converging jet, and (4) the super- 
atmospheric static pressure mound that redirects the merging jet streams in a 
common downstream direction. Comparisons are made between the development 
of the flow, in both the region of jet convergence and the region of combined jet 
flow, and that of the single-jet counterpart which was previously reported. 


1. Introduction 
In an earlier paper by the present authors (Miller & Comings 1957), measure- 
ments of the force-momentum fields in a two-dimensional incompressible free air 
jet were reported. Particular emphasis was placed on the role of the mean static 
pressure stresses in the Reynolds equations of motion. It was shown that the 
pressure stresses and stress gradients are of the same order of magnitude as the 
normal turbulent stresses pu’? and pv’? and their gradients, but of opposite sign. 
As consequences of this near cancellation of stresses in the region of turbulence, 
the deceleration of the mean flow along a streamline depends almost entirely on 
the lateral gradient of turbulent shear stress 7, and mean flow accelerations 
lateral to a streamline are quite small, i.e. the streamlines are nearly straight. 
In this paper a jet flow is considered in which static pressure effects are far more 
pronounced in accelerating the mean flow both along and across the streamlines. 
The two-dimensional flow in question is that produced by the subsonic flow of 
air from twin nozzles set in a common end wall, as illustrated in figure 1. All 
solid boundaries, shown in cross-section, extend sufficiently far in the z-direction 
to ensure two-dimensionality at the plane (z = 0) depicted in the figure. In addi- 
tion, the extension of the end wall (the planez = 0) in the y-directions is effectively 
* Now at Research and Engineering Division, Monsanto Chemical Company, Dayton, Ohio. 
t+ Now at University of Delaware, Newark. 
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infinite. The centre-plane (y = 0) is the plane of symmetry for both flow and 
boundaries. 

Practically no information on dual-jet flows is available in the literature. 
Corrsin (1944) investigated the unstable flow from seven parallel slot nozzles in 
a common wall, with emphasis on flow stabilization methods. His nozzles had a 
width of only 0-20 in., but were geometrically similar to the nozzles of this study. 
Temperature, total head, and flow directions were reported. Bollinger (1950) 
measured mean velocity traverses in the combined flow region of twin round jets 
set with their axes at various angles to each other with the flows converging. The 





Centre- plane 





Fully 
~o| developed 


4 7 
—s ~ flow 
\ ~=—s- Converging flow Combined flow _-__ 


FIGURE 1. Nozzle and flow geometry in the plane z = 0. a = 0-500 in., R = 2-617 in. 





momentum flux distribution in the combined flow region differed greatly trom 
the distribution calculated by adding the momentum flux from the separate jets 
at the same point in space. Considerable curvature of the streamlines was 
indicated. Laurence & Benninghoff (1957) have reported extensive mean velocity 
and turbulence measurements in the fields of multiple interfering air jets pro- 
duced by several multiple nozzle configurations. In all cases studied, however, 
the jets were allowed to mix freely with ambient air, and static pressure effects 
were consequently minimized. 

Despite the paucity of previous data, the major features of the dual-jet flow 
may be inferred a priori from the behaviour of simpler turbulent jets and wakes. 
Thus one would expect to find a region of pronounced negative gauge static 
pressure p between the jets near the end wall, due to the pumping effect of 
turbulent momentum transfer from the individual jets to the central region. 


Consequences of the negative pressure are converging curvature of the individual 
jets and non-conservation of 2-momentum. As the two jets are about to merge, 
one would expect static pressure gradients to be largely responsible for redirecting 
their momenta in the x-direction and for reversing that portion of the flow earlier 
entrained on the concave sides of the jets. After the jets have merged, the 
resultant single-jet flow can be expected to resemble a classical single-jet flow, 
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the similarity becoming more marked with downstream distance. A downstream 
transition from large to small static pressure gradients serves to delineate a 
boundary between the converging and combined flow regions as illustrated in 
figure 1. As might be expected, the mean z-velocity (U) profile achieves full 
development a short distance downstream from this boundary. 


2. Equations of motion 

Under the restrictions placed on the dual-jet flow (and its single-jet counter- 
part), the appropriate equations of motion in customary Cartesian notation, 
using Cartesian co-ordinates, are 


ae cU OV 

(Continuity) = of a= (1) 
Ox Oy 

(z-Reynolds)  —(pU*+pu’*+>p) ee (pUV—7T+yC) = 9, (2) 
ah oy 

(y-Reynolds) 2 (pV2+pv'?+p)+ . (pUV —T—-pE) = 0, (3) 
cy Cx 

where 7 = —pu’v’ is the turbulent shear stress and ¢ = dV /éx—¢cU/éy is the 


mean vorticity. The mean x- and y-components of velocity U and V define a vector 
field identically described by two new scalar field variables—Q, the vector 
magnitude, and /, the vector angle measured counter-clockwise from the positive 
x-axis. While these variables, defined in (7), are hereafter called the mean 
velocity and flow direction, they are not generally identical to the mean magni- 
tude and mean direction of the instantaneous velocity vector. 

At any point in the flow, the above equations of motion may be transformed to 
a local co-ordinate system aligned with the flow by a rotation of the (x, y)-axes 
through the angle 4. In the local co-ordinate system, s measures distance in the 
direction of flow and m (rotated y-axis) measures distance perpendicular to s. The 
transformed equations are 


ee 0G 0 
(Continuity) i AW = Q, (4) 
cs om 
0 C 
s-Rey s — (4pQ? + pul? + p)-—-= — pC) = 0, 5 
(s-Reynolds) a, (BPO + pus” + P)— = (T.— HS) (5) 
(m-Reynolds) Kp@? + = (pv2 +p) iis (7,+ ul) = 0. (6) 
on os 


The subscript s on the turbulent stresses indicates that they are the (s, m)- 
components of the turbulent stress tensor; and x = ¢c@/0s is recognized as the 
local streamline curvature. Several of the important transform relations are 


Qcos6=U, Qsind = JV, (7) 
pu? = pu’ cos? 6 — 27 cos @ sin 6 + pu’ sin? 4, (8) 
pv2 = pu’ sin? 6 + 27 cos @ sin 6 + pv’ cos? 8, (9) 

T, = T(cos* @ — sin® 4) + (pu’? — pv”) cos @ sin 8, (10) 
oQ «=60Ve COU 
ee eee (11) 


cm Cx Cy” 
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Equations (4), (5) and (6) apply to the entire flow field when s and m are 
interpreted as distances along the orthogonal curvilinear co-ordinates formed by 
the streamlines and their orthogonals. We now define the stream function yy which 
satisfies (1) and is zero on the plane of symmetry: 


y= | Udy= | Qcos6dy. (12) 
y 
70 


0 


Streamline slope and curvature are now given by the relations 
dy,/dx = tan0, «x =d0,/ds = cO/es. (13) 


The use of curvilinear co-ordinates in analysing the flow in the converging region 
of the dual jet simplifies comparisons with the single-jet flow. In the single jet, 
the mid-nozzle streamline is a straight line coinciding with the z-axis: in the dual 
jet, each mid-nozzle streamline is highly curved until the individual jets have 
merged into a single jet. In comparisons based on the mid-nozzle streamlines, 
s is equivalent to x in the single jet. A new straight-line co-ordinate n, perpendi- 
cular to the mid-nozzle streamline and tangential to the local curvilinear m- 
co-ordinate, is equivalent to y in the single jet. 

In the jet flows under consideration, all flow variables except the stream function 
must approach zero as y approaches +0. This fact, with the conditions of 
symmetry of the variables, leads to the total mass and momentum integrals 
from (1) and (2): 7 

p | Udy = py, (a constant), (14) 


~ 0 


Px 


| (pU*+pu?+p)dy=J (a constant). (15) 
ev 
A convenient definition of the jet width 6 is 
b(x) = | (U/U,)2 dy, (16) 


«0 


where LU’, is the centre-plane (y = 0) value of U at the x-station of the integration. 


3. Apparatus 

But for the removal of a dummy converging wall originally blocking flow from 
one of the twin nozzles, the flow equipment was identical to that described in the 
earlier paper. A low-speed wind tunnel provided a uniform low-turbulent hori- 
zontal flow of room air to the dual nozzles of aspect ratio 40:1 depicted in hori- 
zontal section in figure 1. The nominal nozzle velocity U, was 72 ft./sec, corre- 
sponding to a nozzle Reynolds number aU;/y of 1-78 x 104. Ceiling and floor 
boundaries over and under the jet flow helped maintain the two-dimensional 
character of the flow while permitting free horizontal mixing of the jets with 
ambient air. 


Mean and fluctuating velocities were measured with a constant temperature 
hot-wire anemometer similar to one described by Laurence & Landes (1952). 
The single tungsten wire, of 0-0002 in. diameter and 0-080in. active length, was 
supported vertically in the flow, i.e. parallel to the z-axis. Measurements of mean 
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static pressure were obtained with a static pressure probe of the disk type and 
connecting Prandtl micromanometer. These same two instruments were em- 
ployed in the single-jet study; a more detailed description appears in the earlier 
paper. 

Two instrument probes were used in measuring and indicating local flow 
direction. The first was a rotatable impact tube supported vertically in the flow. 
A side impact opening, of 0-021 in. diameter, was located near the upper sealed 
end of the tube, which was of 0-065in. outside diameter and 2-5in. length. The 
lower end of the tube formed the journal of a fine half-bearing machined in the 
end of the horizontal supporting member; it also carried a 0-5in. diameter co- 
axial sheave for angular control, and provided a pressure connexion via small 
bore rubber tubing with the micromanometer. A remote graduated hand dial 
and spring loaded piano wire loop communicated known angular settings to the 
impact tube. In use, the direction of maximum indicated pressure was taken 
as the average of two angles of equal indicated pressure straddling the peak 
angle. 

The second flow-direction probe consisted of a light aluminum foil vane, of 
1-5in. vertical length and 0-2in. width, attached along its vertical leading edge 
to a wire spindle carried in the conical end bearings of a fork support. No pro- 
vision was made or needed for accurate vane-angle measurement, since this 
probe was primarily used in regions of high streamline curvature where positional 
accuracy is more important than angular accuracy. 

Measuring probes were individually positioned in the (2, y)-plane, midway 
between the ceiling and floor boundaries, with a precision of several thousandths 
of an inch by means of a lead-screw traversing mechanism. A detailed description 
of apparatus and procedures is available (Miller 1957). 


4. Interpretation of measurements 

Some uncertainty is inevitable in the interpretation of *‘ point’ measurements 
obtained in turbulent flows. Townsend (1956) notes that such measurements 
commonly concern those quantities that are easy to measure rather than those 
which have an easily understood significance. Two reasons may be cited: (1) the 
‘easily understood’ quantities are abstractions of a type precluding direct 
measurement by present techniques, and (2) no existing measuring process is so 
perfectly understood and so free of indeterminate influences as to permit an 
assumption-free definitive formulation of its cause-to-effect relation. The 
confidence placed in a given interpretation, therefore, can only be justified by 
incomplete circumstantia! evidence; examples are incomplete theory, comparisons 
between different instruments, and findings of internal or external consistencies. 
A brief discussion of interpretation practices employed in this investigation 
follows. 

Hot-wire measurements 
The single hot-wire and associated constant-temperature anemometer circuits 
produces a voltage e assumed to be uniquely and instantaneously related 
(through steady flow calibration) to the magnitude q of the velocity vector 
16 Fluid Mech. 7 
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perpendicular to the wire. When the fluctuation of q, i.e. q’, is small relative to 
its mean q, then 
so —=3(4(@)|? =, ; 
f®@=¢=@ and e” afte) = gq? = «3, 
{ de ; 

where f represents the functional e vs q relation obtained by calibration. Measure- 
ments of é and the mean square of e’ were converted to Q and the mean square of 
u, by means of these relations. In flow regions of high turbulence intensity, the 
above expressions no longer hold. Ata point where Q is actually zero (for example, 
at the free stagnation point and the vortex centre), 


f(@) #7 = V(u2+r2) > Q, 


df(é) 


and e”2 
de 


, 
| = q? =ustue—-@ + u. 


These approximate equalities recognize the non-linearity of f at low q (a linearizing 
circuit was not employed). Q values are readily corrected in the neighbourhood 
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FicurE 2. Mean velocity contours. Shaded areas are regions of constant velocity within 
the accuracy of hot-wire measurements. The points marked VC and SP are respectively the 
vortex centre and free stagnation point. The mid-nozzle streamline (y%v = 0-48) starting at 
x = 0, y = 3a is shown, as well as the paths of the four lateral traverses whose profiles are 


given in figures 5 to 8. Values of Q/U,. 


of singular points by interpolation from surrounding regions of low turbulence 
intensity, using the location of the point fixed by vane-probe measurements. 
The method is illustrated in the U/U, profile of figure 10. All computations 
involving Q, e.g. the streamline mapping of figure 3, incorporated this method of 
correction. For purposes of discussion, the Q contours of figure 2 were not cor- 
rected. Turning to the fluctuation term w{?, no correction could be applied since 
no independent estimate of v;? was available. Reported values of this variable 
are subject to considerable error within a radius of 4a of the two singular 
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Flow direction 


The impact direction probe gave measures of the direction of maximum impact 
pressure as contrasted to the flow angle @ previously defined. As is demonstrated 
in figure 3, the two quantities may differ by as much as 40° in a turbulent jet flow, 
where strong lateral gradients of velocity and static pressure exist. Most of this 
error was removed by a method of streamline plotting. This method, discussed in 
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FIGURE 3. Mean flow directions and streamlines determined from the impact probes 
(solid arrows) and the vane (broken arrows). Values of y/aU,. 


conjunction with figure 3, depends strongly on the corrected Q-field and is 
relatively insensitive to impact direction measurements. It is justified by the 
excellent agreement found between the path of the mid-nozzle streamline so 
determined and the path established from impact direction measurements alone. 
Vane directions agreed everywhere with the calculated streamline directions, 
within the ability of the eye to judge. 


Static pressure 

Theory is of no assistance in evaluating the reliability of the static pressure 
measurements, since the external and internal aerodynamics of the probe and 
connecting passage are too complex and obscure. Confidence is instead based on 
the internal and external consistency of the single-jet static pressure measure- 
ments, as presented in the earlier paper. There is as yet no indication of signi- 
ficant bias in the measurements, in either the single- or dual-jet flows. Hinze 
(1959) discusses static pressure measurements at some length. 


5. Results—converging region (0 < 2/a < 12) 

Complete tables and graphs of experimental data are available elsewhere 
(Miller 1957). Measurements revealed, as expected, a high degree of flow sym- 
metry about the centre-plane; consequently, data taken only in the positive 


y half-plane are reported. At x/a = 12, the initially individual jets had merged 
16-2 
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and the flow variables resembled those in a single-jet flow. This was consequently 
selected to be the downstream limit of the converging region and the beginning 
of the combined flow region. 
Mean velocity field 

Figure 2 is a contour map of the mean velocity or Q-field. Contour levels were 
arbitrarily selected to give adequate coverage; and points corresponding to a 
single level were located by traversing the probe through the flow. The value of 
each level was then determined from hot-wire traverses carried out after the 
contours were plotted. To facilitate comparisons with other figures, the positions 
of several characteristic lines and points are indicated in the figure. The reference 
velocity U,, previously referred to as the nominal nozzle velocity, is defined as that 
hypothetical nozzle velocity whose kinetic pressure $pU? is equal to the measured 
gauge impact pressure at the centre of the nozzle mouth. This particular definition, 
selected for convenience in hot-wire calibration and data reduction. does not 
require that Q/U, = 1-00 at the centre of the nozzle mouth. 


Flow direction and streamline fields 


Mean flow directions and streamlines appear in figure 3. Solid and broken arrows 
represent flow directions measured with the impact and vane direction probes, 
respectively. The vortex centre and free stagnation point were found using the 
vane probe. At the free stagnation point the vane whipped around in all directions 
at random without favouring any one; at the vortex centre the vane rotated 
rapidly and continuously. Detection of the two points by the vane method was 
found to be remarkably sensitive and reproducible; no two measurements of their 
positions differed by more than 0-132a or 0-066 in. in either co-ordinate. Crossed 
arrows in the vicinity of the free stagnation point reflect the fact that two sym- 
metric pressure maxima were found on rotating the impact direction probe about 
its axis. 

The mid-nozzle streamline y/aU, = 0-48 was established from impact direction 
measurements exclusively. Using the velocity contours of figure 2, an original 
estimate of its path was made. The impact direction probe was then traversed 
along the estimated path with readings of # taken at 4 in. intervals. These readings 
were used in a numerical integration of the first relation of (13) to obtain an 
improved estimate of the path. Three repetitions of this iterative procedure (the 
last with an x-interval of 0-25 in.) gave satisfactory convergence. Analysis of the 
final curve yielded its curvature « from (13) and are length s as functions of x. 
All other streamlines were calculated from (12) using the measured flow directions 
and Q-values from figure 2. Calculated points are indicated by circles in figure 3. 
In the vicinity of the two nozzle lip streamlines near the nozzle, the agreement 
between streamline and measured flow directions is poor. This is a result of the 
gradient error of the impact direction probe discussed earlier. In these regions, 
the streamlines are more accurate indicators of true flow directions. 

Among the features revealed in figure 3 are: (1) one of the two contrarotary 
stable vortices which recycles air (24°, of that issuing from the nozzles) between 


the jets, (2) the stagnation of the streamline from the inner lip of the nozzle at 
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SP, and (3) curvature, spread, and confluence of the jets. Air entrained on the 
outer side of the jet is seen to approach at an angle near 245°. The corresponding 
streamlines, if sufficiently extended, would reveal another large vortex system 
outside the jet in accordance with the requirement of (14). 


Static pressure field 
Mean static pressure contours appear in figure 4. Comparisons with the previous 
contour maps give an immediate qualitative picture of the important role of 
static pressure gradients in accelerating the mean flow, particularly in regard to 





T T T T T T \ a T T T : me 


y/a 

















” ° y - Fi 72 
FicurRE 4. Mean static pressure contours. Values of p/pU; x 10%. —, contour 
interval = 9; ——-, mid-contour. 


lateral acceleration resulting in streamline curvature. Except in the static 
pressure mound between x/a = 6 and x/a = 10, all regions of the flow were at less 
than atmospheric pressure (negative p). The static pressure mound around 
x/a = 7, whose peak is near the free stagnation point, is largely responsible for the 
stagnation and for the reversal of the flow between jets. Downstream from the 
mound a return to the modest static pressure gradients typical of a single free 
jet is observed. 
Lateral traverses 

Further aspects of the converging flow are best revealed in profile diagrams 
derived from probe traverses indicated in figures 2-4. Proceeding downstream 
along the mid-nozzle streamline, the four lateral traverses perpendicular to the 
streamline and intersecting it at s/a = 2, 5, 8, 12-4 correspond to the profiles of 
figures 5-8 respectively. Traverse measurements included mean velocity Q, 
turbulent s-stress pu?, and mean static pressure Dp. 

The mean velocity profiles of figures 5-7 are fairly symmetrical about n = 0 
and similar in shape (but not of equal peak height) to single-jet profiles taken at 
equal distances from the nozzle. A comparison of the widths of the velocity profiles 
at various downstream distances in table 1 reveals practically identical spreading 
tendencies in the converging jet and the single free jet out to 8 nozzle-slit widths 
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from the nozzle. Beyond 8a the comparison breaks down due to the confluence of 
the converging jets in the dual-jet flow. In table 1 jet width is taken to be the 
lateral distance between the two profile points where the mean velocity is half of 
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FiacureE 5. Lateral traverse at a distance s/a = 2 measured downstream from the nozzle 

along the mid-nozzle streamline. Circles, triangles and squares refer respectively to mean 
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velocity Q, turbulent stress puj?, and mean static pressure p. 0, Q/U,; A, 10u;?/U;; 
KG 12 
O, 4p/pU;. 


Half velocity jet width 


Nozzle width (a) 


x/a or s/a Dual jet Single jet 


2 1-06 1-14 
3 (1-15) 1-20 
4 (1-24) 1:26 
5 1-33 - 
6 (1-46) 1-50 
8 1:96 ——- 
10 2:24 


TABLE 1. Comparison of dual- and single-jet widths near nozzles. 


the profile peak value. Numbers in parentheses were obtained from figure 2 
and are less accurate than the others. Consequences of the nearly equal profile 
spreading are developed later in the discussion of the mid-nozzle streamline 
profiles. 
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FicuRE 7. Lateral traverses at s/a = 8. 0, Q/U,; A, 10u/2/U?; 0, 4p/pU?. 
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Turning to the turbulent stress profiles of figures 5-8, several interesting features 
are observed in the downstream development of the two-peaked curve when 
compared to its single-jet counterpart. (In the single-jet flow, of course, the peaks 
are symmetrical about the x-axis.) Starting near s/a = 5, the inner stress peak 
(negative n) begins a period of rapid growth relative to the outer peak (positive n), 
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Ficure 8. Lateral traverses at s/a = 12-4. 0, Q/U,; A, 10u2/U?; 0, 4p/pU?. 


whereas the outer peak spreads more rapidly and develops a hump clearly 
evident at s/a = 8. Downstream from s/a = 8, in the region where the two jets 
intermingle, the inner stress peak decays quite rapidly; at s/a = 12-4 no trace of 
an inner peak remains. The limited turbulence measurements reported here 
preclude a coherent discussion of the reasons for the observed behaviour of the 
stress profile; it can only be stated that moderate differences exist in the pro- 
duction, transport, and dissipation of turbulence energy between the inner and 
outer shear regions of the converging jet and either shear region of the single jet. 

Static pressure profiles of figures 5-7 bear no resemblance to corresponding 
single-jet profiles. In the discussion of figure 4, mention was made of the role of 
static pressure in curving the streamlines. To explore this further, we rearrange 
(6) in the form ~ ~ - 


D C ( 
4 


C 9 = 
kp? + <— = = (1, + wl) — =— (pv,”), (17) 
cm ma cm 


and hypothesize that the right-hand side is zero along the mid-nozzle streamline 


of the converging jet. Note that in the single jet along the centre-plane each term 
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of (17) is exactly zero, due to the symmetry of the flow. To test the hypothesis, 
values of the two left-hand terms were calculated at various points along the 
mid-nozzle streamline. Resulting values in dimensionless form are listed in 
table 2. These results show that the lateral static pressure gradient at the mid- 
nozzle streamline was at least the dominant force causing streamline curvature. 
That better agreement was not obtained is ascribed more to the diversity of 
measurement errors involved (velocity, flow direction, static pressure) than to the 
influence of any term on the right-hand side of (17). 


s/a KaQ?/U? O(p pu?) €(m/a) 
2 — 0-064 0-055 
3:5 — 0-076 0-090 
5 — 0-076 0-081 
6:5 —0-011 0-011 
6-63 0-000 0-000 
8 0-088 — 0-097 
9-5 0-049 — 0-023 


TABLE 2. Lateral accelerations and static pressure gradients 


on the mid-nozzie streamline. 


Analyses similar to the above have been carried out in other regions of the 
flow. The general conclusion is that both turbulent stress gradients and static 
pressure gradients contribute to lateral flow acceleration; the latter dominates 
where Q is large and xk is small (as along the mid-nozzle streamline), and the former 
dominates where Q is small and « is large (as in the regions of entrainment). 


Longitudinal traverses 

Profiles of Q, ui2, and p, obtained by traverse measurements along the mid- 
nozzle streamline and along the centre-plane, appear in figures 9 and 10, re- 
spectively. Considering the former first, the effect of the static pressure gradient 
along the streamline on the mean velocity is evident. In particular, the static 
pressure peak near 2/a = 8 is largely responsible for the dip in the mean velocity 
decay curve. The turbulent stress profile contains a peak, also near x/a = 8, which 
is sharper and 45 °%, higher than the corresponding peak in the single-jet flow. 

Neglecting the viscous stress term, the Reynolds equation (5) respective to the 
s-direction is ‘ 


OT, 


(4pQ? + pu? + p)— = 0). (18) 


Cs om 
The quantity in parentheses was computed at various stations along the mid- 
nozzle streamline for both the converging and single jets. Results are compared 
in table 3. Good agreement indicates that the distribution of the turbulent shear 
stress 7, about the mid-nozzle streamline was practically identical in the two flows 
despite differences in the mean velocity levels. Equal spreading of the two jet 
flows, discussed in conjunction with table 1, further establishes the similarity of 


the turbulent shear stress fields in the cores of the jets. 
Turning to the centre-plane profiles of figure 10, note the (x, y)-co-ordinates and 
components are used because of the coincidence of the traverse with the x-axis. 
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Plotting U instead of Q reveals more clearly the reversal of the mean flow and the 
nature of the previously mentioned error in mean velocity measurement near the 
free stagnation point. The data of this figure were analysed in a manner similar 
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FiGcurE 9. Longitudinal tranverses along the mid-nozzle streamline. 
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FIGURE 10. Longitudinal traverses along the central plane. 
O, U/U,; A, 10u’*/U?; oO, 4p/pU?. 
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to that of the preceding figure. Equation (18) when written for the centre-plane 


streamline is z 


= (4eU 24 pu’? +p) — = 


y 


= 0. (19) 


a oD ) 


Again the quantity in parentheses was computed at a number of stations along the 
centre-plane. A plot of these data against x yielded a smoothed curve remarkably 
symmetrical about 2/a = 6-57 (very near the free stagnation point) and quite well 
represented by a Gaussian error integral curve as follows: 


w 
(4pU2 + pu’? + p)/pU? = 0-337 | D(w’) dw’ — 0-059, (20) 
—o 
, l rs ee, " 
where O(w’) = —>-\ exp (— 3”), w= 0-658(a/a — 6-57). 
\ (27) 
pQ? + 2pui? + 2p 
pU; 
s/a, x/a Dual jet Single jet 
4 1-00 1-01 
6 0-86 0-87 
8 0:69 0-70 
10 0-55 0-57 
12 0-49 0-49 


TABLE 3. Comparison of mid-nozzle streamline decay. 


Differentiating (20) with respect to 2/a gives the turbulent shear stress gradient, 


ce [(T ; 
——— | —.,] = 0-089 exp [ — 0-216(2/a — 6-57 21 
Ayia (at? pl (x/a— 6-57)"], (21) 
which is a Gaussian error curve centred at x/a = 6-57 with a peak value of 0-089. 


These relations hold in the range 0 < w/a < 11. In view of flow symmetry, 7 itself 
is zero along the centre-plane. 

From the foregoing it is apparent that the turbulent shear stress gradient tends 
to accelerate the mean flow in the positive x-direction at all points along the 
centre-plane out to 2/a = 11; the maximum tendency occurs very near the free 
stagnation point. Alternately reinforcing and opposing this tendency is the effect 
due to the static pressure gradient along the x-axis. In the range 4:5 < a/a < 7-2, 
the static pressure gradient tends to accelerate the flow in the negative x-direction 
against the influence of the turbulent shear stress. Flow stagnation occurs at the 
point (SP) where the two forces are equal and opposite. The x-gradient of the 
normal turbulent stress pu’? in (19) was generally negligible in comparison 
to the static pressure and turbulent shear stress forces. 


6. Results—combined jet region (x/a > 12) 

On merging, the converging jets of the dual-jet flow lose their individual 
identities and behave as a single jet. The analysis of the flow in the combined jet 
region is therefore identical to that of the previously reported single-jet flow 
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(Miller & Comings 1957). Indeed, the combined jet flow is so nearly similar to the 
single-jet flow that it will suffice to describe the similarities and differences without 
repeating the more detailed development of the earlier paper. 


Mean velocity and jet width 
A fully developed turbulent jet is characterized by (1) linear spread with 
downstream distance, (2) power-law decay of the centre-plane mean velocity with 
downstream distance, and (3) similarity of the lateral mean velocity profile. These 
three properties are expressed in the relations 


bla = c(x/a—x,/a), (22) 
U?/U2 = K(xz/a—2,/a)-¢, (23) 
U/U,=f(y), where 4 =y/0. (24) 


In these equations c, x), AK and C are experimental constants, b is the jet width 
(a function of x) defined in (16), and U, is the centre-plane mean velocity at a given 
value of x. Certain restrictions are imposed on the values of the experimental 
constants by the condition of momentum conservation expressed in (15). These 
restrictions, listed in the last lines of table 4 together with values of the various 
experimentally determined quantities, are necessarily approximate since momen- 
tum alone is not strictly conserved, even in the single jet. 


Combined Single 
Factor dual jet jet 
Region of fully 
developed flow 
zla > Ze ij 
c 0-072 0:0723 
xXo/a — 5-50 — 1-572 
K 11-1 6°74 
CG 1:03 1-028 
t(”) exp (— }77*) exp (— }77?) 
Restrictions 
CG 1-0 1-0 
cK = J/pU; = 0-86 0-48 


TABLE 4. Comparison of mean velocity properties in the fully developed jet. 


The data of table 4 with equations (22), (23) and (24) give a complete description 
of the mean velocity fields of the two flows after the mean velocity profiles 
achieved similarity. The station where this occurred, listed on the first line, was 
considerably farther from the nozzle in the dual-jet case because of the extended 
region of convergence of the individual jets. Rates of spread and decay of the 
profiles were identical in the two flows as indicated by the equality of c and C in 
each case. Also identical were the dimensionless analytic functions f(7) best 
fitting the experimental data. The only marked differences found were in the 
position of the ‘ virtual origin’ x) and in the value of the velocity level constant K. 


All comparisons thus far made can be summarized in a single statement as 
follows: if the two nozzles of the dual-jet flow were replaced by a single nozzle of 
the same size as either, located on the centre-plane at z/a = — 3-93 and having 
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a nozzle mouth velocity 1-29 times that of the original dual nozzles, then the 
mean velocity field at x/a > 23 would be indistinguishable in all respects from 
the original field. 

Dimensionless mean velocity profiles appear in figure 11. It is noted that the 
data have not been corrected to remove the slight influence of the lateral velocity 
V on the hot-wire in the region 2-5 < 9 < 4:0. The smooth curve through the 
points obtained at 2/a = 20, 30, 40 is identical to the best curve through the 
uncorrected data of the single jet. Points of the 2/a = 12 profile do not conform to 
the curve since the flow had not reached full development at that station. 
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FicurE 11. Mean velocity profiles for z/a> 12. O, 12; A, 20; 7, 30 and 6, 40. 
The solid curve is identical with that obtained for a single jet. 


Turbulent x-stress and static pressure profiles 
As should be expected, the downstream development of the turbulent 2-stress 
and static pressure profiles resembles that of the single-jet flow. The turbulent 
stress profiles of figure 12 exhibit ‘partial self-preservation’ at z/a = 20, 30 and 40, 
as evidenced by a coincidence of the curves at large 7 and an apparent tendency 
to converge on a single fully self-preserving curve as 2 approaches infinity. 
A convergence of the static pressure profiles on a self-preserving form is apparent 
in figure 13. With the exception of the turbulent stress profile at x/a = 12, all of 
the curves are typical in shape and height to those obtained in a single-jet flow. 
This, however, is as far as the similarities go. No simple equivalence was found in 
the distance histories of the profiles. This is shown in the data of table 5, which 
lists the approximate downstream station of the single-jet flow at which the given 
dual-jet profile was most nearly duplicated. These data reveal a near-normal 
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development of the static pressure field, accompanied by an initially retarded 
but rapid development of the «-stress field. Finally, at x/a = 40, the 2-stress 
profile is ‘in step’ with the static pressure profile; both are nearly identical to the 
corresponding single-jet profiles at x/a = 35. 
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FicurE 12. Turbulent «-stress profiles. z/a > 12: 0, 12; A, 20; 7, 30; O, 40. 
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FIGURE 13. Static pressure profiles. x/a > 12: 0, 12; A, 20; 7, 30; O, 40. 


Thus so far no mention has been made of the nature of the turbulent y-stress 
pv” and shear stress 7 fields. In the earlier paper it was shown that in the fully 
developed single jet (1) dimensionless y-stress profiles are nearly identical to 
coincident static pressure profiles with the sign inverted, and (2) the shear stress 
profile reaches full self-preservation at the same station as the mean velocity profile 
does. We also note that the turbulent x- and y-stresses are identical to twice the 
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a- and y-components of the energy of turbulence. These facts taken all together 
give some insight into the precedence of approach to self-preservation (or dynamic 
equilibrium) of the turbulent motion and energy in the combined jet region. 
First to reach full self-preservation is the turbulent shear stress at x/a = 23. The 
turbulent y-stress or y-energy distribution quickly assumes a characteristic 
single-jet form and approximate growth rate downstream from the region of 
convergence. Lagging behind but rapidly catching up to single-jet proportions is 
the x-stress or x-energy. It is assumed, but without supporting evidence, that the 
dual-jet flow downstream from x/a = 40 would have been found to be completely 
indistinguishable from a single-jet flow originating in a suitably placed single 
nozzle. 


Equivalent single-jet profile 


Dual-jet station, 2/a 
profile po 
station, Turbulent Static 
x/a x-stress pressure 
12 None 13 
20 8 20 
30 22 32 
40 35 34 


TABLE 5. Profile equivalence. 


7. Conclusions 


Based on the experimental data presented and on the results of a previous 
investigation of a single free jet, a reasonably clear description of the mean and 
turbulent motion in a two-dimensional dual-jet flow is possible. In the region of 
jet convergence near the nozzles, the structure of the flow in the core of either 
jet is quite similar to that of a single free jet provided due account is taken of the 
mean flow accelerations attributable to a radically different static pressure field. 
A region of sub-atmospheric static pressure between the converging jets accounts 
for their convergence and for the reversal of a considerable fraction of the total 
flow at the centre-plane against a strong downstream force of turbulent shear. 
Associated with the flow reversal is the appearance of a pair of symmetrical 
contrarotary stable vortices and a free stagnation point on the centre-plane. 
Where the two jets merge, a super-atmospheric static pressure mound redirects 
their momenta into a single resultant jet symmetrical about the centre-plane. 
The flow structure of the combined jet exhibits all characteristics of a single-jet 
flow except for an altered evolution of the turbulence toward self-preservation. 
This is ascribed to differences in upstream conditions. 


For valuable assistance in the course of this work, the authors are indebted to 
Professor Sydney Goldstein of Harvard University and Mr James C. Laurence 
of the N.A.S.A. Lewis Research Center. Financial support was generously 
provided by The Barrett Division, Allied Chemical and Dye Corporation, and by 
the Purdue Research Foundation. 
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The boundary layer near the stagnation point in 
hypersonic flow past a sphere 


By T. K. HERRING 


Department of Mathematics, University of Manchester 
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When a sphere or other bluff body travels at supersonic speeds, a shock wave is 
formed close to the front surface. With increase of speed the air behind the shock 
is further compressed, and the shock wave moves closer to the surface. This paper 
considers the case where the region close to the stagnation point between the 
shock and the sphere can be taken to be a steady laminar boundary layer. 

The approximate solution of the equations of motion follows closely the classical 
work of Homann (1936), ideas similar to those of Lighthill (1957) being used to 
apply it to the problem in hand. It consists mainly in reducing the equations to 
ordinary differential form by assuming forms of the flow variables which satisfy 
the boundary conditions, notably at the shock wave. In addition, several trans- 
formations are employed in order to simplify the equations and to increase the 
range of solutions, and also to facilitate the use of the ‘Mercury’ electronic com- 
puter in solving them. 

The results give an insight into some aspects of hypersonic flows. Included in 
this paper are a selection of temperature and transverse velocity profiles across 
the boundary layer and several graphs relating such quantities as the shock 
stand-off distance and the skin-friction coefficient with Reynolds number. The last 
two mentioned are the most interesting. The first set gives the surprising result 
that the shock stand-off distance increases with increase in Reynolds number, 
whereas it is known that the boundary-layer thickness decreases. The second set 
of graphs shows that the skin-friction coefficient is inversely proportional to a 
decreasing power of the Reynolds number when it is lower than order 10°, but 
the indication is that it tends to the expected constant power of } when the 

teynolds number is of order 104. 


1. Introduction 

In this paper we investigate an interesting case where the region between the 
shock and the sphere can be taken to be a steady laminar boundary layer. A high- 
speed flow through a medium of very low density would seem to satisfy this 
condition. For, at low enough density, the Reynolds number can be so low (say, 
about 104 or less) that the boundary layer will be laminar, but of a thickness 
comparable to the stand-off distance of the shock. 

Now, in the atmosphere at heights of 60-100km, the densities range from 
10-6 to 10-75 g/cm’, and the temperatures vary from 180° to 360° K (Whipple 
1943). The Reynolds number of a sphere of diameter 25 cm, travelling at 4km/sec 
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at an altitude of 80km, is about 5 x 10%. This region of the atmosphere thus pro- 
vides a practical setting for the problem under investigation. 

We shall assume that the surface of the body has not yet been heated signi- 
ficantly, and thus that it is a good approximation to take the sphere and free 
stream temperatures to be equal. 

By confining ourselves to the region close to the stagnation point, the shock 
wave can be taken, to a fair approximation, as spherical and concentric with the 
sphere. This is a result which has been verified under certain conditions both 
experimentally and theoretically. Oliver (1956) showed in experiments with a 
sphere and other bluff bodies at a Mach number of 5-8 that this was so, and he 
included in his paper excellent schlieren photographs of the phenomenon. 
A number of authors, including Lighthill (1957), Freeman (1956) and Chester 
(1956), have put forward theoretical treatments. Chester considered a perfect 
gas with constant specific heats (as we shall), and, by a method of successive 
approximations in which he assumed the shape of the shock and found the shape 
of the body, he was able to verify the above results. 

Now, close to the stagnation point the speeds are subsonic, and all com- 
pressibility effects are due to thermal variations in the fluid properties. The 
present work can therefore be compared with work on low-speed laminar 
boundary layers with varying temperatures. For example, Brown & Donoughe 
(1951) considered the problem of a two-dimensional steady laminar boundary 
layer on a straight porous wall, and they found solutions to the equations of 
motion when the wall temperature was varied. They were able to obtain a negative 
displacement thickness by cooling the wall, a result similar to one in this paper 
where the stand-off distance of the shock is found to increase with Reynolds 
number. 

The geometry of the problem we are considering is close to that of Homann 
(1936), who considered the steady incompressible fiow in a boundary layer, along 
a surface of revolution, close to the stagnation point. He neglected terms of the 
order of the square of the transverse co-ordinate measured along the surface, 
and found that the form of the flow variables suggested by the axial symmetry 
and the constant thickness of the boundary layer, which is a correct result in this 
approximation, satisfied the boundary conditions and the equations of motion. 

In the paper by Lighthill (1957), an exact solution to the inviscid incompressible 
high-speed flow past a sphere was obtained. Lighthill showed that all the equa- 
tions of motion and boundary conditions could be satisfied assuming a spherical 
shock concentric with the sphere. 

For the present problem we combine the approaches of Homann and Lighthill, 
using the forms of the flow variables suggested by the boundary conditions behind 
the shock wave. The use of the strong-shock conditions is found to be necessary 
to make the density dependent upon the normal co-ordinate yet independent of 
the transverse co-ordinate to the first order. With this simplification, the assumed 
forms of the flow variables satisfy the boundary conditions accurately. The 
inviscid terms of the equations of motion are satisfied accurately, and the viscous 
terms are satisfied to the first order in the transverse co-ordinate. A note on the 
approximation involved in using the strong-shock relations is included in §3. 
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We will assume for simplicity that the fluid is a perfect gas and that the specific 
heats at constant pressure and volume of the gas mixture are constant throughout 
the region of interest. In a real flow, however, where there are very large tem- 
perature differences in the boundary layer, the variations in the ratio of the 
specific heats can be as much as 0-3. This is attributable to ionization and dis- 
sociation effects as well as to the large range of temperatures. For this reason we 
shall quote results for different ratios of the specific heats. 

The first section of the paper is devoted to the formulation and simplification 
of the equations of motion, and several transformations are introduced to 
facilitate subsequent integrations. The second section concerns the derivation of 
the results and the discussion of interesting physical conclusions. 


2. Basic equations 


The equation of state of a perfect gas is 
p = pRT, (1) 


where F is the gas constant per gramme, p the pressure, p the density, and 7’ the 
temperature. To simplify the problem further, we will take the viscosity ~ to be 
proportional to 7’, and the Prandtl number o = yc,/k to be a constant = 0-72. 
Here c,, and ¢, are the specific heats at constant pressure and volume (c,/c, = y), 
and k is the thermal conductivity, which must also be proportional to 7’. 

We shall use spherical polar co-ordinates (7, 4, A) with origin at the centre of the 
sphere and the axis 6 = 0 pointing upstream, and write the respective velocity 
components as ¥,, Vs, V,, With v, = 0 in the present case of axial symmetry. The 
radius of the sphere is taken to be a, and that of the shock to be b, and suffixes 0, 1 
will refer respectively to the uniform conditions ahead of the shock, and those 
immediately behind. Also, we will denote the enthalpy by J and let M, = Vi/eo, ep, 
and V, be the values of the Mach number, speed of sound and axial velocity ahead 
of the shock. It is important to note that throughout this paper we are considering 
only the region where @ is small (shaded in figure 1). 

The strong-shock relations are 


Pp = y+! (2) 
Mh t= 
y—1.. aa 
Mn = — 2 Wo008 9, tm = sind, (3) 
Da 
Pi _ _*Y  f2c0s70, (4) 
Po yrttl 


ail I, = T, fy _ 24W- 1) M2 


= - - 5 
Ih TM fo (y+1)° ) 


These shock boundary conditions suggest that it may be useful to investigate 
whether an approximate solution with 
v.=v(r)cos#, v, = u(r)sin§#, (6) 


p=plr), p= p(r)+P(r)sin?6, (7) 
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may exist. The more complicated form of p is particularly suited to a boundary- 
layer type of theory, which is well known to be sensitive to the effect of pressure 
gradient. The assumptions (6) and (7) correspond to those of Homann (1936) and 
Lighthill (1957). 

The equations of motion in spherical polar co-ordinates are reduced to the 
axisymmetric boundary-layer form by dimensional considerations and are listed 
below. The continuity equation becomes 


= (pv,) + 


—— = sin @ ) = Q, bad 
cr asin@c6é (sin (pv) (3) 


and the momentum equation for the 6-direction becomes 


CUg Upp, lop c 4) (9) 
Cr}? 7 


1) eel A ne ean i 
P\’r or a cb, acO or 


Sphere 


—_——> 
nee ene ‘J 
6=0 > 





Uniform region 





ie 
Bg 
FIGURE 1 


whilst the energy equation reduces to 


el vy el Vy Cp le/ ol at 
Rial ne a | C ; 
P(r, aco} acé 5 ar! 5.) +H( 5: (10) 
where ee 2 (11) 


In accordance with the usual boundary-layer theory, the pressure gradient 
across the layer is neglected, so that, in (9), ¢p/e@ is found from its value im- 
mediately behind the shock, which is assumed to be within the boundary layer. 


Thus (4) gives On 4y 
Py Ss / Mm sin @ cos OD: 
Al) (y+1) 


2 2 


and by replacing 1/2 by V2/c? and then c? by ypo/po, this becomes 


1. Sah. 
20 ~~ (y+ip V2, sin. (12) 


On substituting for 0p/c0 = Cp,/c in (9), we get 


| Up | ve 0) 4(iy-1) V2, 0 [ *e) 
LL +— an | = sinf+ —{u—}]. 13 
Pra, a co (y+1)* a — +5, (¢ or (19) 


To compare this work with that of Homann, expressed in the boundary-layer 
co-ordinates (x,y), we simply replace a@ and asin@ by x, and r—a by y. On 
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substituting for the flow variables the forms (6) and (7), we are able to reduce the 
equations of motion to ordinary differential form. Equation (13) becomes 
(du wu*\ A(y—-—1)V2p, d / 2) ; 
v—+—}— ~ — — - Ss ?0= 0: 14 
lo( dr a ) (y+1)?a dr (x dr}\*" (14) 
and we shall also use (10) for 9 = 0, namely 


(15) 


Pr. ld (5*) — 


dr odr\' dr 
where henceforth J(r), 7'(r), 4(r) will signify values for 0 = 0. 
To satisfy the continuity equation, we now choose a Stokes stream function y, 


such that ous ays 
asin Opv, = —p, 50 and sin Opv, = p; an} (16) 


and on substituting (6), we see that iy takes the form 


y = V(r)sin? 9, (17) 
and (16) reduces to 
2 », TY 
v= fiw and «=O. (18) 
ap p dr 


Transformation of the equations 


The effect of compressibility is made less marked by choosing a Howarth variable 


n=c | Pag, (19) 
Ja Pi 
C being a constant. We put also 
V(r) = Bf(n). (20) 
where £ is another constant, and 
I(r) T(r) _ Pi (21) 


g() = = ap ed 2 
’ I, T’ p(r) 
Equations (18) now become 


u(r) = 


2 p 

a p(7) 
where the prime signifies differentiation with respect to 7; and substituting (22) 
in equation (14), we get 


Bf(y), u(r) = BCf'(y), (22) 


aqy2 2072 y—1) V2 
_9 B ( ff" + B _f'2 ~ 4(y 2 H 09— d af") BCy, =) 
a a (y+1)? a dy \'4P, 
Now, since “p = const., if we choose 
/ ry 1 
C= ( ro.) and B= (Vv,a)}, (23) 
ya ali 
the above equation reduces to 
: - ee 
"4.9 (a 24. ~g = 0, (24) 
ee he dae 


and the energy equation (18) becomes 
g" +2ofg’ = 0, (25) 


where the Prandtl number o0 = 0-72 for air. 
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The equations (24) and (25) represent a fifth-order system, but in addition the 
stand-off distance of the shock is unknown. We require therefore six boundary 
conditions to solve the equations. These are listed below. 

First at the sphere 7 = 0, the no-slip conditions give 


v,=01,=0, or f(0) =f*'(0) = 9; (26) 
and since insignificant heating is supposed to have occurred, 7'(0) = 7); whence 
from (5) sa 1, _ yt 1)2 f as 

T, 2y(y-l) IR 


The remaining three conditions are found at the shock, 7 = 7, say. From the 
shock relations (2) and (3). and using (22), we have 
Viiy-1 2 
w.=- 0) pur Bf (n,), 
(y+1) a 


where B = (V,v,a)}: and substituting for v, this becomes 


4V2a*y(y—1)? M?2\3 

B= | + 1 ) 
2V,a V, 
where Ry = Vou Po and M, =o 
Ho Co 


are the Reynolds number and Mach number. Hence we find 


sony = 7 (Pee. (28) 
Also u, = Vo BOS"(g,) = BS'(s,), 
whence tS '(a,) = 1; (29) 
and 7'(y,) = T,, whence 9(9,) = 1. (30) 


Using a simplified coding technique for a machine such as the Mercury 
electronic computer, the equations can be integrated from the sphere along an 
outward normal. Five boundary conditions at the sphere are needed to do this, 
and so two additional conditions on f"(0) and g’(0) are prescribed, approximate 
values being guessed from a crude polynomial solution. 

As the problem stands, the shock-wave position is determined when f’(7) and 
g(7) are simultaneously equal to 1. It would therefore be necessary to try many 
runs with varying values of f”(0), g’(0) and g(0) before this could be accomplished. 


Further transformations 


To get over the above difficulty, we first put 


G = : = q. (31) 


so that (24) and (25) become 


f+ 2" $246 = 0, (32) 
and G” + 2afG’ = 0, (33) 
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whereas the boundary conditions on G are 


G(0) = a (34) 
| fea 

4(y—1) ‘ 

G v0 


A shock-wave position is now determined when f’(7,) = 1, and represents the 
solution for that value of y which satisfies (35). Unfortunately, y, M, and R, are 
now determined uniquely by the boundary conditions (28), (34) and (35), and 
thus only one shock-wave flow is found for each set of initial conditions. This 
transformation then offers little improvement in the solution. Our final trans- 
formations, however, enable us to determine several flows for each integration. 
First we choose a new variable 
a = Ay, (36) 
where A is a constant, and then by putting 

St (9) = Ax(a), Gy) = Aty(a), (37) 

we are able to retain (32) and (33) in the same simple form: 
xv" + 2xx"—2'2+y = 0, (38) 
y" + 2oxy' = 0, (39) 


where the primes refer in this case to differentiation with respect to a. 

The boundary conditions, however, now include an arbitrary constant A, and 
y, M, and R, are not fixed by the boundary conditions, but depend on the value 
A takes. 

We now have the conditions at a = 0 


2(0) = 2'(0) = 0, (40) 
» 

and y(0) = —_.. (41 
K\ yM2A4 ) 

At ~ = a, say, where a, = Ay, we have 
‘ Pe Ry y+ ; 9 
“(a4) = (a9 [12 (42) 

‘ 1 
(a) = 5» (43) 
y= 4) 

1 = : 44 
anc Yy(%,) (y+ 1)? r4 ( ) 


It should be noted that real solutions of (43) and (44) exist provided 2’(a,) > 0 
and y(a,) < 0-5, and that the smaller of the two solutions of (44) for y is chosen on 
physical grounds. 

The two most important advantages gained by these transformations are: 
(1) The boundary conditions can now be solved uniquely for y, M, and R, at each 
integral step (subject to the above conditions). (2) There is no loss in generality 
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if we fix y(0) and thus reduce the initial conditions to be varied from three to 
two. A convenient value useful in the polynomial solution is seen to be 0-02. 
Another advantage is the added usefulness of the crude polynomial series 
method for finding suitable values of x”(0) and y’(0), since no longer does this type 
of solution have to apply to specific values of y, M, and R,. In fact it is only 
necessary to guarantee that 2’(ax) and y() increase steadily with « to assure 
satisfactory solutions of (43) and (44) and thus of (41) and (42) also. 


3. Presentation of the results 

The tables 
The results obtained from the Mercury electronic computer by integrating 
equations (38) and (39) for 60 evenly spaced values of x”(0) and y’(0) were printed 
in the form of the tables 1 to 4, and others which are included in the author’s 
thesis (Herring 1958). 


2 m4 4 5 6 ie 8 9 10 

y(p)ydp x(a) y(x) u’(a) x”(a)/y(a) y(a) y'(a)/y(a) 

Jo 

x 0 0-000 0-000 0-000 30-000 0-020 12-500 

a= 02 0-009 0-001 0-119 8-445 0-070 3°568 

a= 0-6 0-057 0-018 0-347 3°221 0-169 1-436 

a= 1-0 0-143 0-075 0-550 1-776 0-262 0-831 

0-264 0-184 0-717 1-093 0-340 0-505 

0-412 0-340 0-849 0-730 0-397 0-290 

hoe 0-579 0°527 0-953 0-545 0-433 0-147 

0-756 0-728 1-041 0-464 0-451 0-063 
a 0-937 0-938 1-122 0-435 0-458 0-021 Y My R, 
0-1158 1-121 1-156 1-201 0-427 0-460 0-006 1-497 9-82 4,86) 
0-1020 1-305 1-386 1-280 0-425 0-460 0-001 1-407 10-79 7,70 
0-0908 1-489 1-629 1-358 0-425 0-461 0-000 1-343 11-72 11,60 
0-0816 1-673 1-886 1-436 0-425 0-461 0-000 1-294 12-63 16,81 
0-0738 1-858 2-158 1-514 0-425 0-461 0-000 1:255 13-52 23,50 
0-0672 2-042 2-444 1:593 0425 0-461 0-000 1-225 14:39 32,26 
0-0615 2-226 2-745 1-671 0-425 0-461 0-000 1-200 15°25 43,11 


and 6 are indicated by *. Rows in this table only 


0-4 except for the first step. 


In each table the columns are numbered from 1 to 10. The columns 2-7 are 
functions of x, chosen to be directly proportional to the velocity and temperature 
variables, and columns 1, 8, 9 and 10 give the shock stand-off distances and the 
values of y, M, and R,. The rows are printed at intervals of « = 0-2 except in 


are printed at intervals of a = 


table 1 where intervals of « = 0-4 are used for brevity, and those numbered 
represent the positions of the shock waves. It must be remembered that each 
table gives a picture of the flow across several completely independent boundary 
layers. This is clear from the relationships connecting the columns, which are 
functions of x, with those expressed in terms of r. The constants connecting them, 
which are different for each shock, are inserted in square brackets, in the following 
expressions. 

Column | gives the ratio of the stand-off distances of the shocks to the sphere 


radius, namely b—a 


= —  y(B) ap, 
0 


a Sar(ax,) x(a) J 


(45) 
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4,86 
7,70 
11,60 
16,81 
23,5 
32,26 


43,17 


only 





om whee 





b-—a 
a 
0-0945 
0-0821 
00725 
0-0648 
0-0584 
0-0529 
0-0482 


a 
0-1037 
0-0921 
0-0826 
0-0747 
0-0679 
0-0620 


a 
0-1013 
0-0844 
0-0724 
0-0632 
0:0558 


9 


a 


y(f) dp 


J 0 


0-000 
0-009 
0-028 
0-057 
0-095 
0-143 
0-198 
0-261 
0-330 
0-404 
0-481 
0-560 
0-641 
0-723 
0-806 


TABLE 


9 


+ 


0 
0-000 
0-017 
0-060 


TABLE 3. The shock data used in figures 4 and 8 


” 


. 


0-000 
0-017 
0-060 
0-128 
0-220 
0-334 
0-467 
0-614 
0-773 
0-938 
1-108 


y(B)dB x(x) y(x) 


y(p)dB x(a) y(a) 
0 


ax’(a) 


0-000 
0-159 
0-316 
0-468 
0-613 
0-751 
0-880 
1-001 
1-114 
1-221 
1-323 
1-421 
1-518 
1-614 
1-709 


ax’(a) 


0-000 
0-239 
0-472 
0-693 
0-900 
1-089 
1-260 
1-413 
1-553 
1-682 


4 
x’ (a) 


0-000 
0-279 
0-552 
0-813 
1-061 
1-291 
1-503 
1-700 
1-886 
2-063 
2-235 


0 


xv” (a)/y(a) 


40-000 
11-303 


6°451 
4-417 
3°291 
2-580 

“099 


-766 


v 


ax" (x)/y(a) 


0-000 
7-891 
4-091 
2-660 
1-900 
1-433 
1-127 
0-925 
0-794 
0-714 
0-668 
0-644 
0-633 
0-628 
0-627 


Jv 


xv” (x)/y(x) 


70-000 


9-226 
4-812 
3°164 
2-301 


ee ee 
bo 
w 


0-997 
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6 
y(a) 


0-020 
0-070 
0-120 
0-168 
0-215 
0-259 
0-298 
0-331 
0-358 
0-378 
0-393 
0-403 
0-408 
0-412 
0-413 


2. The shock data used in figures 3 and 7 


6 
y(a) 


0-200 
0-150 
0-279 
0-404 
0-522 
0-628 
0-717 
0-789 
0-841 
0-876 
0-898 
0-910 
0-916 
0-919 
0-920 


6 
Y(&) 


0-020 
0-150 
0-279 
0-403 
0-519 
0-621 
0-705 
0-769 
0-814 
0-842 
0-857 


od | 


1-423 
1-351 
1-296 
1-253 
1-218 
1-190 
1-166 


are indicated by *. 


“474 
“396 
-336 
-289 
+251 


990 


1 daa 


vw 
é 


“480 
“375 
+304 
-251 


-209 


TABLE 4. The shock data used in figures 5 and 9 are indicated by *. 
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whereas column 2 gives 
ee 
y(P) dp. 


70 





This can be converted to (r —a)/a upon multiplying by [(y + 1)/82(a,) v’(x,)], but 
this is more quickly found from (45) as G a— v/| ‘pyap| by dividing the 
#0 


corresponding values in columns | and 2. Note that the suffix 1 denotes the values 
immediately behind the shock and these are given in numbered rows. 

The next five columns, 3-7, are connected to the velocities and temperature as 
follows. Column 3 gives (a) y(«), where the ratio of the radial component of the 
velocity behind the shock to that ahead is — v(r) cos 0/V, cos 6, where 


u(r) y-1 
= | - | a |) ferate) x(x) y(a). (46) 
0 ie 
The transverse component is quickly found from column 4. which gives 2’(«). 
| : 5 
since (r) 1 
u(? ‘ ™ 
") Ee | (2). (47) 
0 1 


To obtain [du(r)/dr]/Kj/a. column 5 must be multiplied by [8x(,)/(y+1)]. 
The temperature ratio is 


T(r) y(%) 

= + = [50] y(a), 48 
T, = (0) [50] y(«) (48) 
where (a) is given in column 6, and finally column 7 gives y’(«)/y(a), where the 
constant needed to give [d7(r)/dr]/[T,/a] is 


8: a rf 40 
Fee ed E ~ w(a)2'4)|. (49) 


Certain information can be found by direct reading of the tables. For instance, 
the edges of the velocity and temperature boundary layers are approached when 
columns 5, 6 and 7, which are related to du/dr, T and dT'/dr, tend to a constant 
value. This is very clearly seen in, for example, table 1 (see also figure 2), and in 
this, as in a few other cases, the shocks are well outside the boundary-layer 
section. The assumption of a constant pressure throughout the region behind the 
shock (12) will be less accurate in these cases. 

In the table mentioned above the thickness of the velocity and temperature 
boundary layers are seen to be of the same order, and this is due to the nearness 
of the Prandtl number (0-72) to unity. 


The graphs 


Three temperature and two transverse velocity profiles have been drawn for 
tables 1-4. Figures 2-5 give the temperature and figures 6-9 the velocity profiles 
for each of these tables, taken in order. The magnitude of the temperatures 
reached depends directly on the shock relations and a comparison of figures 2 
and 4 gives some idea of the range considered. 
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FicurE 2. Temperature profiles plotted from data given in table 1. N.B. In figures 2-9 
a shock wave is denoted by a stroke | and approximate boundary-layer edges by a cross x ; 
the number above is the curve number. Curve 1: y = 1-200, M, = 15-25, R, = 43,170. 
Curve 2: y = 1-294, M, = 12:63, R, = 16,810. Curve 3:y = 1-407, M, = 10-79, 
RK, = T7700. 

FIGURE 3. Temperature profiles plotted from data given in table 2. In curves 1 of figures 3-5 
the shock-wave positions and the boundary-layer edges suggested by the profiles have been 
marked. In actual fact of course, the boundary layer will end at the shock wave, since 
ahead of it the temperature assumes the constant value 7'5, and so outside the shock the 
curve is shown as a dotted line. Curve 1: y = 1-190, My) = 14-79, Ry = 6340. Curve 2: 
y = 1-296, M, = 11-62, R, = 1900. Curve 3: y = 1-423, M, = 9-34, Ry = 660. 
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FicurE 4. Temperature profiles plotted from data given in table 3. Curve 1: y = 1-220, 
M, = 20-55, Ry, = 22,070. Curve 2: y = 1-289, M, = 17-96, R, = 11,270. Curve 3: 
y = 1-396, M, = 15-27, R, = 5030. 
Figure 5. Temperature profiles plotted from data given in table 4. Curve 1: y = 1-209, 
M, = 20-32, R, = 9750. Curve 2: y = 1-304, M, = 16-51, R, = 3480. Curve 3: 
y = 1:375, M, = 14-50, R, = 1850. 
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Naturally in this problem we are most interested in, for example, the shock 
stand-off distance, or the boundary-layer thickness as a function of the Reynolds 
number #, or Mach number ,. To find the variation with #,, we must maintain 
y and M, at constant values. 

A comparison of curves | in figures 2 and 3 and also curves | in figures 4 and 5, 
where in each pair y and M, are approximately constant, gives an indication of 
this behaviour. The approximate boundary-layer edges have been marked on 
these graphs with a cross and in these examples as R, increases the gradients get 
larger and thus the boundary-layer thickness decreases. However, at the same 
time it is surprising that the shock stand-off distance (marked with a dash) 
increases with #,. This last result is borne out in a more detailed investigation 
later. 
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Figure 6. Transverse velocity profile plotted from data given in table 1. Curve 1: 
y = 1-200, M, = 16°25, Ry = 43,170. Curve 3: y = 1-407, M, = 10°79, K, = 7700. 
FicurE 7. Transverse velocity profile plotted from data given in table 2. Curve 1: 
y = 1-190, M, = 14-79, R, = 6340. Curve 2: y = 1-423, M, = 9-34, Ry = 660. 
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FicurE 8. Transverse velocity profile plotted from data given in table 3. Curve I: 
y = 1-220, M, = 20-55, R, = 22,070. Curve 3: y = 1-396, M, = 15-27, R, = 5030. 
FicurEe 9. Transverse velocity profile plotted from data given in table 4. Curve 1: 


y = 1-209, M, = 20-32, R, = 9750. Curve 3: y = 1-375, M, = 14:50, R, = 1850. 
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In figures 6-9 showing the transverse velocity profiles, the edge of the boundary 
layer is defined by a region of approximately constant shear flow. In other words, 
in a region where viscosity has little effect the vorticity shed by a shock wave of 
constant curvature is constant. The outer parts of these graphs agree quite well 
with the exact inviscid incompressible solutions obtained by Lighthill (1957) for 
a gas with varying specific heats. The effect of viscosity in reducing the velocity 
quickly to zero is very marked. 

The remaining graphs, figures 10-16, show the variation of certain quantities at 
constant values of y and M,, and provide a very interesting insight into their 
behaviour in a hypersonic flow. 

Interpolation on a large scale is an unfortunate necessity in this method, where 
y, M, and R, are not predetermined, and to maintain y and M, constant, a three- 
point linear scheme between two tables was used. 
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FIGURE 10 FIGURE 1] 
FiGuRE 10. The Reynolds number as a function of the ratio stand-off distance of the shock 
to sphere radius, when y = 1-4. Curve 1: M, = 15. Curve 2: M, = 12. Curve 3: M, = 9. 
FicurE 11. The Reynolds number as a function of the ratio stand-off distance of shock to 
sphere radius, when y = 1-3. Curve 1: M, = 18. Curve 2: M, = 15. Curve 3: M, = 12. 


Figures 10-12 show the variation of R, with the ratio of stand-off distance to 
sphere radius for a number of constant values of y and M). The interesting result, 
hinted at previously, that the shock stand-off distance increases with Reynolds 
number, is here verified, and the asymptotes to the graph (which are independent 
of Mach number) are drawn for each value of y. Can this mean that we have 
negative displacement thicknesses? A comparison of figure 2 with 3, and figure 4 
with 5 gives an indication that the boundary-layer thickness decreases with 
increase in Reynolds number and an examination of the product pw in these 
cases shows that it decreases from a maximum at the shock to a minimum at the 


6 
sphere. In other words, the displacement thickness given by [ (u,p,— U(r) p(r)) dr 
“a 


is positive. The shock distance, however, is not governed by the boundary-layer 





270 T. K. Herring 


thickness, for in table 1 the shocks are outside, and in tables 2—4 they are inside the 
boundary layer, and thus it is probable that this idea of displacement thickness is 
erroneous for the region behind the shock. If the shock position depended more 
upon distribution of mass than of mass flow, i.e. depended more upon p and pu, 
“b 
then a ‘mass defect’ given by | (p,—p(r))dr would certainly be negative, and 
7a 
the shock distance would increase with increase in the Reynolds number. 

A rough interpolation within figures 10-12 to maintain constant y, and between 
them to maintain constant J, shows that the shock distance decreases with 
increase in Mach number and with decrease in y, respectively. A simple analysis 
of the boundary conditions on the density at the shock and the sphere, and the 
knowledge that the density decreases continuously outwards from the sphere 
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FicuRE 12. The Reynolds number as a function of the ratio stand-off of shock to sphere 
radius, when y = 1:2. Curve 1: My = 18. Curve 2: M, = 15. Curve 3: M, = 12. 


provides a simple explanation of these results. The strong shock condition on the 
density (2) is independent of Mach number, whilst the density at the sphere, 
which, since 7'(0) = 7, is dependent only upon the pressure given by (4), in- 
creases with Mach number. The region therefore is further compressed and the 
shock wave moves closer to the body as the Mach number increases. A similar 
consideration of (2) and (4) shows that as y increases, the density just behind the 
shock decreases, whilst at the sphere it remains effectively constant. The shock in 
this case, therefore, moves away from the sphere as y increases. 

A note here about the error caused by assuming the strong shock conditions 
on the density will be of value. The oblique shock relation for the density is 


y+1 1 


M1 2-3 
— = > 50 
po ae ae (8) 


where terms O(/?) have been neglected. The strong shock relations are only valid 
if 4.M2(y—1) > 1 which means, for example when y = 1-4, that 1/2? must be much 
greater than 5 for these relations to be used. The lowest Mach number considered 
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at y = 1-4 is 9, and substituting in (50) gives p,/p) = 0-94(y + 1)/(y—1). It is 
clear from this that the strong shock conditions tend to overestimate the density, 
and so the correction to the graphs involves a small shift to the right to give a 
compensating increase in the shock stand-off distance. This correction will of 
course decrease as the Mach number increases. 

Figures 13 and 14, where the ratio of (r—a) at 7’ = 0-987, to (r—a) at T = 7), 
i.e. (b—a), is plotted against log,) po, give an indication of the variation of the 
boundary-layer thickness with Reynolds number. The decrease in the ratio with 
increase in #, shows clearly that the temperature gradients are increasing, and 
this supports the idea that the boundary-layer width is decreasing with increase 
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FicuRE 13. Boundary-layer thickness as a function of shock stand-off distance, as indi- 
cated by the ratio [(r—-4@)7.o987,]/(0—4@), for y = 1-3 and M, = 18. 

FiGuRE 14. Boundary-layer thickness as a function of the shock stand-off distance as 
indicated by the ratio U(r — 4) 70.987, ) (6—a), for y = 1-2 and M, = 15. 


in Ry. At low Reynolds number, where the profile approximates to a straight line, 
the ratio approaches 0-98, and it decreases through a point of inflexion before 
tending to a further constant value for large values of Ro. 

Finally figures 15 and 16 show the variation of skin friction with the Reynolds 
number, where the skin-friction coefficient is defined as 
Moldu(r)/dr],., _ 4002"(0) 


C= 5 ; 
, bPo VE (yiy+)D}EM, 


and it is noticeable that, with y and M, held constant, c, is solely dependent upon 
the initial value of x”(~). Once more the influence of the shock wave on the flow 
can be clearly seen, for, instead of c,R} having a constant value as suggested by 
boundary-layer theory, it is seen to decrease with increase in Ry. In fact, when 
R, is of order 10%, c, behaves more like Rj *. However, the graphs seem to indicate 
that c,R} does approach a constant value at higher Reynolds numbers, and this is 
not surprising since the tabulated results indicate that when A, is of order 104 the 
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shock waves lie well outside the boundary layer, and thus one would expect the 
normal boundary approximations to apply. 

The approximate solutions which have been found for the problem have only 
been explained in part, but further clarification should be possible as more 
becomes known about flows at very high speeds. 
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FicurE 15. The skin-friction coefficient as a function of the Reynolds number, for y = 1:3 
and M, = 18. 
FiaurE 16. The skin-friction coefficient as a function of the Reynolds number, for y = 1-2 


and M, = 15. 


I would like to thank Prof. M. J. Lighthill and Mr C. R. Illingworth for their 
help and encouragement at all times, and also Mr R. Kerr for his advice on 
programming for the Mercury computer. 


REFERENCES 
Brown, W. B. & DonovuaueE, P. L. 1951 Nat. Adv. Comm. Aero. Tech., Wash., Note 2479. 
CuHEesTER, W. 1956 J. Fluid Mech. 1, 353, 490. 
FREEMAN, N.C. 1956 J. Fluid Mech. 1, 366. 
Herrina, T. K. 1958 M.Sc. Thesis, Manchester University. 
Homann, F. 1936 Z. Angew. Math. Mech. 16, 153. 
LIGHTHILL, M. J. 1957 J. Fluid Mech. 2, 1. 
Outver, R. E. 1956 J. Aero. Sci. 23, 177. 
WHIPPLE, F. L. 1943 Rev. Mod. Phys. 15, 246. 











he 


ly 
re 


Wn 





A presentation of cnoidal wave theory for 
practical application 


By R. L. WIEGEL 


Department of Engineering, University of California, Berkeley 
(Received 2 May 1959) 


Cnoidal wave theory is appropriate to periodic waves progressing in water whose 
depth is less than about one-tenth the wavelength. The leading results of existing 
theories are modified and given in a more practical form, and the graphs necessary 
to their use by engineers are presented. As well as results for the wave celerity and 
shape, expressions and graphs for the water particle velocity and local accelera- 
tion fields are given. A few comparisons between theory and laboratory measure- 
ments are included. 


1. Introduction 

Mathematical arguments show that the theory of surface waves, commonly 
known as Stokes’s waves, is most useful (i.e. valid without unreasonable restric- 
tion as to the wave amplitude) when the depth to wavelength ratio d/L is greater 
than about 1/8 or 1/10 (Keulegan 1950; De 1955). For shallower water the theory 
for a wave type known as cnoidal appears to be more satisfactory. The formula 
for the wave profile involves the Jacobian elliptic function cn w; hence the term 
‘enoidal’, analogous to ‘sinusoidal’. Keulegan has pointed out that the validity 
of this theory rests on the assumption that the square of the inclination of the 
water surface is small in comparison with unity. (On the other hand, for Stokes’s 
waves to be valid in shallow water, the wave amplitude is required to be exces- 
sively small, thus making the theory unrealistic.) The theory of cnoidal waves has 
not been developed in the literature to the state where it can be used extensively 
by engineers; and the object of the present paper is to fill this need. 

Korteweg & de Vries (1895) initiated the theory of cnoidal waves. They showed 
that it accounts for a general class of long waves of permanent type and finite 
amplitude; one limiting case of the theory gives the solitary wave, while another 
limiting case gives sinusoidal waves as accounted for by linearized wave theory. 
Keulegan & Patterson (1940) have studied the cnoidal wave on somewhat 
different lines. Again, Keller (1948) treated the problem using the general non- 
linear shallow-water theory and obtained formulae which are similar to those of 
Korteweg & de Vries. Littman (1957) has proved the existence of permanent 
periodic waves of the type in question. The approximate region of validity of the 
existence proof of cnoidal waves by Littman is shown in figure 1. Benjamin & 
Lighthill (1954) have advanced the theory considerably with regard to the 
formation of bores and hydraulic jumps, and Iwasa (1955) has also considered it. 

18 Fluid Mech. 7 
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A comparison between the various alternative developments of the theory will not 
be attempted here; rather, we select those results which appear to be most useful 
and proceed from them to derive data for practical application. 


Sub- super- 
V{[H/d]K(k)}=0-50 


~— wn itic: , 
e : / 





H/d 











0 1-0 
C?/gd 


Figure 1. Approximate regions of validity of existence proofs by Struik for Stokes’s 
waves and by Littman for cnoidal waves (after Littman 1957). 


2. Presentation of the theory 

Kortweg & de Vries, Keulegan & Patterson, and Keller use different symbols; 
however, the critical formulae obtained by them are essentially the same. The 
wavelength is given by 


L 4 oe Y; —$ 
q= gh) (2241-7) (1) 


where d is the still water depth, A(X) is the complete elliptic integral of the first 
kind with modulus & (it should be noted that A(k) is sometimes denoted by 
F(k)). y, is the vertical distance from the ocean bottom to the wave trough and 
L and k are defined by the following two equations: 
(Y, ld) - (y/ d) 
2L+1 —(y, d)’ 


(2041 — 4) E(k) = (21 L+2-%_ my) K(k), (3) 


d 


2 = 


where y, is the distance from the ocean bottom to the wave crest and E(k) is the 
complete elliptic integral of the second kind with modulus k. The following 


inequalities must also hold: 


- Y; ° 
2h+1>:°: 2 and O0<f? <1. (4) 


In using tables of elliptic integrals and functions, the reader is warned that they 
are often tabulated as functions of the parameter m, where m = k?. Equation (2) 


eiateadiaieinaee (22 “a a _ Ye a) —(y d) - H (5) 
d k? ke 
Substituting this into equation (1) and squaring gives 
L?H ‘ 
— = 18(kK(k)}*. (6) 


di 
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FIGURE 2, Curves obtained from equations (6), (10) and (116), showing relationships 
between L?H/d° and k*, K(k), and y,/H —d/H = y,/H —d/H +1. 
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Fiecure 3. Curves showing relationships k? vs L?H/d*, and k? vs T(g/d)t and H/d. 
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Both Keller (1948) and Littman (1957) also obtain this relationship as their basic 
approximate solution. L*H/d* is plotted as a function of ‘? in figures 2 and 3. If 
the wavelength, wave height and water depth are known, then the many formulae 
of the cnoidal wave theory can be used as they are expressed in terms of various 
functions of the square of the modulus k. The terminology of the elliptic functions 
and integrals as used herein are as used by Milne-Thompson (1950). 

The wavelength is 


(sy 


3H ) an) (7) 


Equation 3 can be rearranged to give 


. ’ (1—[y,/d]) ,, 
E(k) — K(k) = — ; K(k). 8 
: (2 +1—[y/d]) 6) 


Substituting (1) into the above equation gives 


y. 16d? ; 
T= gpa (K(R)[K() — B(R)}} + 1; (9) 
or, multiplying by d/H, we get 
y. ad 16d 
: k) (K(k) — B(i)]}. (10) 


HH ~ 3H 
Next, y, can be obtained from the relationship 


y y. H _ 16d? H 


—-—= (k)}}+1——, 

aa ae = g72 (K(h) )[A(k) —E(k)]} + 1 7 (lla) 

y, ad 16d* | — } 

or HH * oe apt ( k) (K(k) — E(k)}}. (115) 


These equations have been plotted in figures 2 and 3. The relationships among k?, 
K(k) and E(k) have been tabulated over the range £2 = 1-10~* by Kaplan (1946, 
1948) and partially by Hayashi (1930, 1933) and Airey (1935). In order to extend 
these functions to the range needed for the study of waves (4? = 1—-10~*°) the 
following equations were used (Jahnke & Emde 1945): 


K(k) = A+ }(A—1)k'2+3,(A—3) k'4+ 25 (A-30) ko +..., (12a) 


256 
E(k) = 14+3(A—4)k’2 4+ 3,(A—43) kh’ + JS (A-S) ko 4+ .... (125) 


where k’ = ,/(1 —k?) and A is the natural logarithm of 4/k’. 


The wave profile is 


y, = ¥,+H en? 2x) (7 -~ 7) . X ; (13) 


where cn is the Jacobian elliptic function associated with cosine. The function en 
is singly periodic ~ ided k is real number and 0 < k < 1. The period becomes 
infinite when k = | (in which case we have the sehen wave). While the period of 
the en function is 4A (xk), the period of the cn? function is 2K(/). The expression 
en? [2K(k) (z/L—t/T), k] is plotted in figures 4 and 5 as a function of 2/L, t/T, 
with k? as parameter. Values of the cn function are available over a limited range 
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of k? (Milne-Thompson 1950; Spenceley & Spenceley 1947; Schuler & Gabelein 
1955). Values for 0-95 < k? < 1 were calculated to three decimal places using the 
following series (Milne-Thompson 1950): 
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FIGURE 4. Surface profiles of cnoidal waves, obtained from equation (13). 
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FicurE 5. Surface profiles of cnoidal waves, obtained from equation (13). 
where wu is the incomplete elliptic integral of the first kind. x/L is the same as 
u/2K(k) for the en? function. uw has been used rather than the more commonly 
accepted symbol uw to avoid confusion with the horizontal component of water- 
particle velocity. In figure 6 the theoretical surface profile is compared with some 
measurements made by Taylor (1955). It can be seen that the cnoidal theory 
predicts the wave profile very well. 

It is interesting to note that when the modulus k£ is zero, 


en(u|k) =en(w|0)=cosu and K(k) = 37; 
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hence 4K(k) = 27, and the wave profile is given by the trignometric functions. 
When &= 1, en(u|1)=sechu and we have the hyperbolic function with 
K(k) = «; hence the period becomes infinite and we then have the solitary wave. 
It can be seen that when k is reduced from | to 60-9999 the period 4K is reduced 
from infinity to about 77, whereas the further reduction of k to 0 reduces 4K to 
only 27. 
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FicuRE 6. Comparison between measured and theoretical wave profiles. 


The wave velocity (using Stokes’s second definition of wave velocity, which is 
the velocity of the propagation of the wave-form when the horizontal momentum 
of the liquid has been reduced to zero by the addition of a uniform motion) is 





ro HLL E(k)Y ' 
t; = ( 4 +— ann j 

(gd) i 1 BA mE) | (15a) 

—_ 16d3 "2):) = l 
i a 

16d? 4. (1 E(k) ' 

C = (¢ } T 2(k _ . BL 

(94) ! sre KM) (5 a vais 


Equation (15a) has been plotted in figure 7. 
For one limiting case (the solitary wave), k? approaches unity, (4) is unity and 
K(k) is infinity; hence 


H 
C= (gd)i (1 +5, . (15c) 
This approximation to 
c= (oa (1+2)\" (15d) 
r ok. Se 


is higher by a maximum of only 2° even for the case of the solitary wave of 
maximum steepness (H/d = 0-78). 

For the other limiting case (the linear theory, where k? > 0), E(k)/ K(k) > unity 
and ; 


972/J2 
C = (gd) (1 a 7) (15e) 
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which is an approximation to 





























































































































































































































































































































































































































































































































































































































































































































































































































4n°d2\)\3 
Owe elt 15 
\9 302 | ( f) 
Now, the linear theory gives (Lamb 1945, p. 366) 
‘ 4 
gL 27d\* . 
Gin — tanh ——} . 15 
| 2Q7r L | (159) 
eetoh the 
2 ini 
‘= Suman +0-50 i A 
eu annus iii HS 0-45 a 
f: eS ke iat tt 
iiipn penance sis 
wr A mm RELL 
aoaa nasa ana 
es oe 
penetest ition ars 
kW Ba RSE | 
Ht 
rs: BaaaE j 
i jaaee 
nif +tpiit + pee beets ui 
i i } 
i if T 
pauatt 
Ni an 
pauen 
isnt 
sigan - 
aieauie ui 
in 
1000 10,000 
LH /d3 
FIGURE 7. Curves obtained from equation (15a), showing relationship between the 
cnoidal wave velocity (Korteweg & de Vries) and L?H/d*. 
But the first two terms in the expansion of tanh (27d/L) are 
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which is in agreement with equation (15 f). 

The equation for wave celerity as given by Keulegan & Patterson (1940), and 
Littman (1957), which is the velocity of the wave crest with respect to fixed 
co-ordinates, can be written 


ti cite 1 2-35) ; 
Ot = gil +5 |-1+35(2 3 Kb) I: (16) 


This equation has been plotted in figure 8 as a function of L?H/d°, and in figure 9 
as a function of T(g/d)}. 
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FicurE 8. Curves obtained from equation (16), showing relationship between the 
enoidal wave velocity (Keulegan & Patterson) and L?H/d*. 
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FicurE 10. Comparisons between theoretical and measured wave velocities. 
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In figure 10 are shown comparisons between the wave velocity as obtained 
from equation (16) and some measurements of waves in the laboratory. 
As cnoidal waves are periodic and of permanent form, | 
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Figure 12. Relationships among T(g/d)!, L?H/d? and H/d (Keulegan & Patterson). 


and the wave period is given by 
4 2 I\ 4 IK (1 
m(g\? _ (16d\? kkK(k) 7 
ra) = (sx Hi Bp ial 
+ ae \2~ K(k) 





using the velocity as given in equation (15a). 
From this equation k? can be determined as a function of 7T'(g/d)? and H/d, and 
from this L?H /d* can be determined. This has been done and plotted in figure 11. 
Using the velocity as given in equation (16) results in 


r(2)' " (sir) = kK(k) . (18) 
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which has been plotted in figure 3 as a function of k? and H/d, and in figures 12 and 
13 as a function of L?H/d? and H/d. The results in figure 13 show that a wave of 
a given period and height, and for a given water depth, can have two possible 
lengths. The physical significance of this is not apparent. 

The pressure at any distance y above the bottom has been shown by Keller 
(1948) to be, to the second approximation, 


P = PIYs—Y); (19) 


where p is density and y, is given by (13). It is rather surprising that this simple 
hydrostatic expression applies. 

The horizontal and vertical components of water-particle velocities at any 
point x, y within the fluid can be obtained from the following equations given 
by Keulegan & Patterson (1940): 


u fh Pe d )< a 
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nfs ” ; | amis i 
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In order to use the equations for water-particle velocities the necessary 
numbers can be obtained from figures 2, 3, 4 and 5, as 


u|k = 2K(k)a/L, (24) 
sn? (u|k) = 1—en? (wz | hk), (25) 
dn? (u | k) = 1—k[1—en? (uw | &)). (26) 


A few comparisons of theory with laboratory measurements are shown in 


figure 14. In considering the vertical velocity it should be noted that the curve 



































t/T 
oa 0-25 050. 
i | | | 2 
— | ~ 
sg le 
—— S 
s we. = 
& = = oo. 
a aoe 
3 yj oon = 
5 -20 ee i | -20 8 
oO -40 4 4 4 4 i 1 i 4-10 - 
a 0 0-5 1-0 15 20 25 30 
t (sec) 
0-4 1 Z : 
§ 04 = Crest L ; 
03}——+ + 03/SWL/ “f 
\ Trough |— F % 

g 02 7 is a 

rl \ / rs 
o1;}-—1+ 01-47-14 

olL_i a us 
0 02 1-2 0 O02 04 06 O08 


lv| (ft./see) 


FIGURE 14. Comparisons of horizontal components of water-particle velocity and accelera- 
tion with cnoidal wave theory. a: 0, Experimental points (Elliott 1953); - - - - - , Stokes’s 
second-order theory ; , Cnoidal theory. b to d: e,o, Experimental points (Morison & 
Crooke 1953); ———, Linear theory; - - - - - , Stokes’s second-order theory; , Cnoidal 
theory. In } and c, points e,o give |w| under wave troughs and crests respectively; in d, 
points e, o give |v| under SWL leading wave crests and following wave crests respectively. 








H (ft.) T (sec.) d (ft.) L (ft.) 
a 0-483 3°20 2-0 24-7 
b 0-105 1-62 0-292 5-10 
c 0-120 1-27 0-292 3-71 
d 0-137 2-09 0-292 6-58 


plotted for the cnoidal theory is for the phase of maximum vertical velocity and 
this occurs prior to the time the wave profile goes through the still water level. 


The work presented in this paper was done under Grant G-4630 from the 
National Science Foundation. The author wishes to express his appreciation to 
Mrs M. M. Turner for her help in the computations and preparation of the 
illustrations. 
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Optimum power generation from a moving plasma 
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The possibility exists of directly using the plasma, resulting from a controlled 
fusion reaction, to generate electricity by electromagnetic induction. Two special 
cases of a more general problem are considered here: (1) the extraction of optimum 
power from the steady one-dimensional flow of an incompressible inviscid plasma 
across a uniform transverse magnetic field in an externally loaded channel of 
arbitrarily varying cross-section, and (2) the extraction of optimum power from 
the steady one-dimensional flow of a compressible inviscid plasma across a uni- 
form transverse magnetic field in a channel of uniform cross-section. In each case, 
the magnitude of the required external loading at optimum power operation is 
determined as a function of the parameters which characterize the hydro- 
magnetic interaction. Also determined are the magnitudes of the terminal voltage, 
power, fluid mechanical to electrical conversion efficiency, and the variation of 
the fluid dynamical variables along the channel at optimum power. 


1. Introduction 


With controlled thermonuclear fusion a possibility for the future (Post 1956), 
the question may be raised as to the best, most efficient way of utilizing the energy 
liberated for electric power generation. The conventional steam-turbine method 
of generating electricity, which involves moving parts of heavy cumbersome 
machinery, may not in this case be the most practical and efficient way of electric 
generation. A possible method of electric generation which does not involve 
moving parts is generation by electromagnetic induction using the highly con- 
ducting plasma of the reaction products as the working fluid. The inductive action 
may be described as follows. A plasma moving perpendicular to a magnetic field, 
both flow and field being horizontal, has induced in it a vertical electric field. If 
the flow is in a channel with top and bottom walls conductors, and these walls are 
connected by an external load, then current will flow through the plasma and 
externalload. Further, itis contemplated that by utilizing some of the mechanical 
energy of the exhausted ionized gases in a plasma pinch engine (Kunen 1958), 
electric power may be generated to be used either as a prime or auxiliary power 
unit for the engine system. 

The dynamical state of the plasma, characterized by the state variables uo, po, 
p, and 7, exiting from a thermonuclear reactor or exhaust of a pinch engine will 
probably be fixed by the operating power level of the reactor. Also, the channel 
entrance width and length may be fixed by space and other limitations. We then 
pose the following question. For a channel of given entrance cross-section and 
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length, what must be the shape of the channel, distribution of applied magnetic- 
field strength and magnitude of applied external loading so that maximum 
power may be extracted by the external load? The derivation of the conditions 
for maximum power transfer from conventional generators is developed in 
elementary texts, but for a plasma generator of the type consid>red here the 
situation is much more difficult to analyse because of the complex hydromagnetic 
interaction between a variable magnetic field and a compressible electrically 
conducting fluid continuum. 











Ficure 1. Diagram of plasma generator. 


To simplify the analysis, consider the two-dimensional channel shown in 
figure 1. The channel is of unit breadth out of the paper. The plasma moves in the 
positive x-direction as shown, entering the channel of fixed opening y) with fixed 
values Uo, Po, Pp and 7, of the state variables. A magnetic field, normal to the 
paper and directed into it, is assumed to be an unknown function of x only. The 
following simplifying assumptions are made: 

(a) The flow is one-dimensional, i.e. the fluid dynamical state variables vary in 
the x-direction but not over the cross-section. 

(b) The magnetic Reynolds number (see § 5) is small; that is, any effects on the 
fluid flow of secondary magnetic fields resulting from the induced current dis- 
tribution are negligible, because either the secondary magnetic fields are too 
small or are in the wrong direction to produce appreciable effects. 

(c) The induced currents in the x-direction are small compared to those 
induced in the y-direction, and are neglected. 

(d 

(e) The intrinsic properties ¢ and y of the plasma are taken to be constant. 

(f) The fluid is assumed inviscid and non-heat conducting, so that the only 


The plasma behaves like a perfect gas. 


~— 


dissipation is electrical. 
(7) The plasma is assumed to be electrically neutral so that no space-charge 


sheaths are developed near the conducting walls. 
(h) The motion of the plasma through the channel is smooth, i.e. shock free. 
(t) The channel walls are perfectly conducting. 
The extremum problem, as outlined above, was formulated in Neuringer (1958), 
using the techniques of the calculus of variations to obtain the appropriate system 
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of differential equations and boundary conditions. These equations are each of 
first order and highly non-linear, and they require machine calculation for their 
solution. In the special cases where the differential constraints, i.e. the fluid flow 
equations, are integrable, analytical solutions may be obtained, however, without 
resorting to the mathematical apparatus of the calculus of variations. The two 
special cases considered in this paper are: 

(1) The flow of an incompressible, inviscid plasma across a uniform transverse 
magnetic field in a channel of fixed length in which the variation in cross-section 
along its length is arbitrarily prescribed. 

(2) The flow of a compressible, inviscid plasma across a uniform transverse 
field in a channel of fixed length and uniform cross-section. 

In each of these special cases, we shall determine the magnitude of the required 
external resistance, at optimum power generation, as a function of the parameters 
which characterize the hydromagnetic interaction. Also, we shall determine the 
resulting terminal voltage, power, fluid mechanical to electrical conversion 
efficiency, and distribution of the fluid dynamical state variables along the channel 
at optimum power. 


2. Notation 


The rational M.K.S. system is used throughout. 


u = velocity k = voltage appearing at the terminals of 
p = pressure the generator 

T = temperature n = efficiency 

p = density m = constant mass flow rate through 

E = electric intensity channel 

3 = magnetic induction c, = speed of Alfvén waves 

H = magnetic intensity P = dimensionless pressure 

J = current density (’ = dimensionless velocity 

I = total current X = dimensionless distance along channel 
r = resistance per unitlength &’ = dimensionless terminal voltage 

R = total resistance M, = Mach number of flow at channel 

ao = plasma conductivity entrance 
ft, = magnetic permeability a = Ily, 

x = distance along channel 6 = magnetohydrodynamic interaction 

y = channel width parameter 

1 = channel length x = parameter (function of y and M,) 

y = ratio of specific heat Pf = parameter (function of y and M,) 

P =. elactvie power () = subscript representing entrance values 


3. The circuit and magnetohydrodynamic equations 

Let r(x) represent the external unit resistance at station 2 along the channel. It 
may be defined as that element of the parallel resolution of the total external 
resistance which would carry the full current density associated with that station. 
Let J(x) represent the induced current in the y direction per unit length of channel 
at stationa. The external voltage drop is then J(x) r(x). The voltage drop due tothe 
19 Fluid Mech. 7 
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internal resistance of the plasma is J(x) y(x)/o. Kirchoff’s law then gives for the 
closed circuit at station x: 


Jy 
Ey = Jr+—, (1) 
where E = uB. (2) 
Assumptions 7 and c imply Jr = k, (3) 


where k, the voltage drop across the external load, is a constant to be determined. 
Using (2) and (3) in (1), we obtain 


- 
J = . {uBy —k}. (4) 
The power generated in the external circuit per unit length is 
P k 
=rJ*=kJ = iiled | 
length \ y) 
and the total power generated is 
el / i) 
P = ke | (ut - = dx. (5) 
0 ¥ 


The problem then is to find the velocity distribution u(x) and the constant k 
which when inserted into (5) maximizes P. 

The one-dimensional compressible fluid flow equations, modified to account for 
the hydromagnetic interaction, are (Resler & Sears 1958): 


Continuity equation , 
PUy = PoUpYo = M, (6) 


where m represents the constant mass flow rate through the channel. 
Momentum equation 


du _P 


—_ 
dx dx D 


=-—-JB= — (ou? — ; (7) 


where JB (with proper sign) represents the Lorentz force per unit volume 
exerted by the magnetic field on the fluid. The last term is obtained using (4) for J. 
Energy equation 

d | 

dx | 


2 
Fe i+ ip | uy! i {puy} + ud By — . y = 0, (8’) 
where the first term on the left-hand side represents the net flux of internal plus 
kinetic energy through the faces of a volume element, the second term represents 
the rate of mechanical work done by the pressure forces, the third term represents 
the rate of doing work by the Lorentz force, and the last term represents the 
Joulean dissipation. (8’) can be written, by means of (4), as 


d | 
dx | 


— 
v0 


Vy k 
L__© 4 dat ia tholuB— | = 0. 
y-lp * "y 1 
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An interesting observation may be made using the energy equation in the 
form (8). Substituting for the integrand in (5) using (8) and performing the 
integration, we obtain, using (6), 


Rae oe ee 
I = (5243 +a] nl Laas tae | (9) 


That is, the power delivered to the external circuit is equal to the difference in 
the total, or stagnation, enthalpy flux at the entrance and exit of the channel. 


4, Case 1: incompressible inviscid plasma flowing across a uniform 
transverse magnetic field in a channel of arbitrarily prescribed cross- 
section 


Let y(x) represent the arbitrarily prescribed variation of channel cross-section 
with length. The fluid density is uniform along the channel, so that the fluid 
velocity, from (6), is given by u(x) = m/pyy. Substituting for the velocity in the 
expression for the power (5), we obtain 


A * rl 
| gee -— -ia) | da . (10) 


For maximum power, we require dP/dk = 0. Carrying out the differentation, 
we obtain 
Bm UpYoB 
x oe (11) 


2Po 2 


Two interesting observations may be made using (11). 

(1) The value of the external or terminal voltage drop at maximum power is 
independent of the channel shape. 

(2) Using (11) in conjunction with the condition that the induced e.m.f. at each 
station along the channel equals uwBy = u,By, = constant, it is seen from Kir- 
choff’s law that the internal voltage drop is equal to the external drop and hence 
a local matching theorem for maximum power is valid; i.e. the applied external 
unit loading, r(x), at station z must be equal to the unit internal resistance of the 
working plasma at that station. It is not too surprising that in the case of incom- 
pressibility local matching is required for optimum power extraction. This is so 
because incompressibility yields the condition that the induced e.m.f. at every 
station along the channel is the same. This, coupled with the fact that the external 
voltage drop at each station must be identical (the channel walls are perfectly 
conducting), means that each station is uncoupled from the others and so can be 
treated as independently isolated. Conceptually, we may visualize the incom- 
pressible plasma generator as one consisting of a continuous distribution of 
independently acting elemental generators each with its own internal resistance. 
The elements (generator plus internal resistance) are then connected in parallel 
and the resulting combination feeding energy to an external load connected to the 
conducting walls. 


19-2 
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The maximum power is, substituting for & from (11) into (10), 


P = uayg to [ “=, (12) 
JoY(*) 

The analysis leading to (12) was based purely on electrical energy considera- 
tions, i.e. Kirchoff’s law and the definition of electrical power in terms of current 
and resistance. It is instructive, therefore, to reconsider the optimum power 
formula from a fluid mechanical energy viewpoint. Before doing so, the pressure 
distribution along the channel, at optimum power operation, is required. Sub- 
stituting the value for & obtained in (11) into the momentum equation (7), we 
obtain the following differential equation for the pressure distribution along the 


channel: 


dp m* dy lobB?m 
dx pyy'dx 2 poy” 


integrating. and using the boundary condition p = py at x = 0, we obtain 


m* {1 1 | ob?m [7 dé 


P= Pye 5 : ‘ 
' “poly Yi! 2Po Jo ¥(S) 


Now the fluid mechanical energy density of an incompressible, inviscid fluid is 
given by the Bernoulli expression (p + $p,u?). The energy flux at any station 2 is 
then m/po(p+4p,u*). Energy conservation requires that the difference in the 
fluid mechanical energy flux at the entrance and exit of the channel goes into 
electrical energy. Since we have already demonstrated that at optimum power 
operation the matching theorem is valid, the electrical energy must consist of two 
equal parts: (1) the electrical energy extracted by the external load, and (2) the 
electrical energy delivered to the fluid and which is dissipated internally in the 
form of heat. Mathematically, we have 

P= = up + Pot |y—[p 7 5Po u* |}. (14) 
~/o 

Substituting the expression (13) for p into (14), we find that P reduces identi- 
ally to (12). 

Incidentally, since the fluid pressure can not be less than zero at the channel 


exit, we have 


m= { | 1) oB?m [' dé | 0 (15) 
— pee ae : is ” 
Po 2/ My \y?(L) yt 2) . oS) 


as the condition which must be satisfied between the fluid mechanical, electro- 
magnetic, and geometrical properties of the system. 

Finally, as an interesting check on the inner consistency of what has gone before, 
we shall derive the maximum power formula using the techniques of lumped 
circuit theory. The external conductance in length dz is dx/r(x). The total conduc- 
tance is then l dx 


nt 


and the total external resistance is 


f de} 
Rext 75 {| . 
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The current flow in the external load is 


rt 
= [ Jayde, 
0 


; da 
and since J = k/r(x), we have r= | =, 
01 (x) 
The total power delivered to the outside circuit is /*Rext. Using the above 
expressions for J and Rext, we obtain 


P = I*Rex, =k? [S. (16) 
0 


Now at optimum power, k = }(uyy)B) and r(x) = y(x)/o (local matching). Sub- 
stituting for k and r(x) into (16), we see that (16) reduces to (12). 

For later comparison with the compressible flow case, we repeat the principal 
result of this section: for an incompressible plasma, the condition for optimum 
power requires that the external load be equal to the internal resistance of the 
plasma. 


5. Case 2: compressible inviscid: plasma moving across uniform trans- 
verse magnetic field in a channel of uniform cross-section 
Analysis 
We treat the optimization problem for the compressible flow case in a different 
manner than the incompressible case. Instead of seeking to maximize the power 
integral (5), we shall avoida complicated integration by considering the maximiza- 
tion of the equivalent expression for the power, namely (9). (9) states that the 
electric power delivered to the external load is a maximum when the total 
enthalpy flux at the channel exit is a minimum. However, before we can mini- 
mize the total exit enthalpy flux, we must first integrate the fluid mechanical flow 
equations (6), (7) and (8), in order to determine the appropriate expressions (to be 
inserted into (9)), for the velocity, pressure and density of the plasma at the 
channel exit. 
Before integrating, it is convenient to introduce the following dimensionless 
variables: 
Ue £. p 
Uy Yo. PoXs 
and the following dimensionless parameters: 


, _ UY%B I Po I Lb? "yp 7 al 
i’ = aa =a; 3 = ype? 
k Yo Pots yMG m Yo 
The meaning of the first three parameters is fairly obvious and needs no further 
consideration. Let us briefly consider the significance of the fourth. Remembering 
that B = n,H and the definition of m, we can write é as 
Byol op? Hl 1, HH? 
6=—* ae = = (1,019) (© A} 
m Yo Poo Poo 
In the next to the last form, 6 represents the ratio of the electrical body force to 
the dynamic or inertial force on the conducting fluid (Resler & Sears 1958). Hence 
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when 4 is large we can expect the applied magnetic field to produce significant 
effects on the fluid motion. 

Consider now the resolved form of 6. The term (,alu,) is defined as the mag- 
netic Reynolds number (based on channel length). It is a measure of the effect of 
the hydromagnetic interaction on the magnetic field, and may be interpreted as 
the ratio of the motion induced magnetic field strength to the applied magnetic 
field strength. The term (yu, H?/p,u2) can be written as 


Sa > 
Pole uz’ 
where ¢,, = (44, H?/po) represents the speed of propagation of certain magneto- 
hydrodynamic waves, called Alfvén waves. Thus ,/(p)ug/,H*) is analogous to 
a Mach number where the usual speed of sonic disturbances is replaced by the 
speed of the Alfvén wavelike disturbances. The dimensionless number, 6, we shall 
henceforth call the magnetohydrodynamic interaction parameter. 


In terms of the dimensionless variables and parameters, equations (7) and (8) 


become dU dP 6é(,, 1 
MNT an coe a ; 
dX *dX" a | v) ie 
d Y 6 / l 
—_ PU+4U? U- = Q. 8 
misc? +3 +a 7) m7 
Eliminating (U — 1/k’) from (7’) using (8’), (7’) becomes 
dU dP et + vial 
ey vig _— LU?) = 0. 
ax*ax—* ax es ‘ened | 
Integrating, using the boundary conditions U = 1 and P = 1/yMj at X = 0, we 
obtain ’B) —(U —(4k’) U? 
p. G+) -W-(P)" (17) 
[ony 1)]k 
1 (% 4. 
where desi + Me’ =? or Mz" 3) 


Substituting for P in (7’) using (17), we obtain, after some differentiations and 
combination, 


dU _ 4 eee | az 
7 | (I-52 ¥ 4) = (1-80) (1-% #0) 


y ee 6 / 
dp KY(a+kp)—(t s ,u)}| + (1- 
york 2° /)J avy 
Integrating (18), using the boundary condition U = 1 at X = 0, we obtain the 
following expression for U’ as an implicit function of X: 


EE ee ee 1—k'U\ 
oe telat a teal adie kp) log( 7) 


y— ye+1)\(y—1 Pp 
D | deat )( + ly =I) k (a +K'p)} log ( 


“i ny Day 
y zy 
/ / 


a ev) (v- D, =0. (18) 


—[y/(y—1)]#’U 
l—[y/(y—-1)]F’ , 


yp spepy Ltt pl | 5 
tn 6+") >" lap epi a 
(19) 


ea \ 3, 
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The total exit enthalpy flux, in terms of the dimensionless variables, is 


Y 
mus f pP@ U(a)+4l xa), 
where the argument a is the value of X at the channel exit. The necessary 
condition that the total exit enthalpy flux be a minimum requires the derivative 
with respect to k’ of the above expression to equal zero, or 


y {dP(a) dU(a)) dU (a) 


(y—1)\ dk’ a (+ Ul) ay = © (20) 


U(a)+P(a) dk’ 


Differentiate (17) with respect to k’ to obtain dP(a)/dk’. Substituting for 
dP(a)/dk’ into (20), rearranging and simplifying, we obtain, 


+B=0, (20’) 


adi _ (y= 1) (a+k'Z)-AY-1I V+ yly+ KV? 

where A= Ay —1)?{1—-[y/(y-1)] kU}? 
pg 27-1) BU + 27°aU?— (+1) US 

(vy 1 {1—[yl(y— 1) POF 


Differentiating (19) with respect to 4’, rearranging and simplifying, we obtain 


14243} 5 44546474849} = 0, (22) 
where 
og ED ot cat ie 
1= PH 2 I= y(y— 1) ka +k A); 
— yk’ {l y+] ). 
= + yk" B 
2 fly DIU) ly By FEA 
he "(a +k’ B)—(y?-1)k} 
2(y— 1) (1—[y/(y— 1) kU)” 
4 = —y(y—1) (2+ 2k’) log ( A ); 
1—[y/(y—1)] k’U 
§ = y(y—1) (a+ 2k’f) )log (~, DI ee 
(l—WU) 


| 2 } (a+ PY eae)? 
(1 Hy? +1) +y%"(a + k/A)} 1 — UO} 
~ (l=[y/(y-V kU) = [y/(y- ky 

_ f(y 1) Rak! + 3k) — My + D}(U 1), 


B= Ul bily— FO) (I= bity— DIF) 
9 _ WAlly- Dk a+ ie Sy + DE ylty— IW - {+ U)-ylty— 1K DY, 
7 —[y/|(y-I) kK UP =[yl(y—1) kK’)? 
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Eliminating dU’ (a)/dk’ from (20’) and (22), we obtain 
1+2+3)}B-—{4+5+6+7+8+9A = 0. (23) 


Equation (23), together with equation (19) (evaluated at X = «), form a system 
of two simultaneous transcendental equations in the two unknowns k’ and U(a) 
as functions of the parameters y, M, and 6. 


é 


Results and discussion 
A specification of the ‘input’ parameters y, M,, 6 and a, and the quantities kh’ and 
U(a) obtained from the simultaneous solution of (23) and (19) is sufficient to 
determine completely all of the electrical and fluid mechanical properties of the 
interaction. In particular, we shall now derive, in terms of these parameters. 
explicit formulas for (1) the ratio of the required external resistance at optimum 
power to the internal resistance of the equivalent incompressible flow, (2) the 
power, and (3) the generator efficiency. 
Eliminating J from (1) and (3), using (2), and solving for 7, we obtain 


Kyo 


~ o(uBy,—k) 


Proceeding similarly as in $4, we sum in parallel the continuous distribution of 


these elemental resistances, and obtain 


a pt _s 
R e kyo | | (wBy, —k) ia 
eT 


ext 
/0 


An element of internal resistance is y,/o. The total resistance, Rint, of an incom- 
pressible plasma moving in a channel of uniform width yp is then y)/ol. Forming 
the resistance ratio, integrating the second term of the above integral, and 


transforming to dimensionless variables, we obtain 
Rext Ce” : fk (¥@_dZ ,.. I 
na | ip: oe | geo a, ie (24) 
Rint a Jo : a Ji dl 
Forming dX /dU as an explicit function of U by differentiating (19), substituting 
into the latter integral and performing the integration, we finally obtain 
Rint 


> 
Rext 


k’ 
= s [10+ 11 + 12] — 1, (25) 
a i) dee! Ce me ee A 
where 10 = oF “3-70 — 1) (a +2'A)) log ( ~~ ): 
Sg 
2k'y Pl i-pyly-D# J 
| yk! (a + kB) — (y? — W/2y _47+F), 
(1 -[y/(y-DIRFU)Q-[yy-DIR) oy 


Bin iy ~1)(a+k'f)- 


12 =(U-1) 





where it is understood that U’ is evaluated at X = a. 
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We emphasize once again that in the above formulas, Rint does not represent 
the internal resistance of the compressible plasma, but the resistance which would 
obtain had the plasma been treated as incompressible. It is clear that only in the 
case of incompressibility, where the induced e.m.f. at each station is the same, can 
we add in parallel the continuous distribution of elemental internal resistances. 
We chose to normalize Rext with respect to the incompressible Rint because the 
latter is a resistance property of the system which is independent of the hydro- 
magnetic interaction and depends only on the conductivity of the fluid and the 
geometry of the system (i.e. the width to length ratio). 

The formula for the optimum power is obtained as follows: 


P = I?Rex, = k?| —— (see § 4). 


A [ dx 
J o7(%) 


Substituting the expression for 7(x), obtained above, into the integrand, and 
transforming to dimensionless variables, we obtain 


p - Puoryo! E [" UdX-1 | 
k’ aJe ° 


The bracketed term is precisely (Fint/ ext); hence 


~ Bruroyol Rint 


P= a , 26 
k _ Text ( ) 
We define the efficiency or effectiveness, 7, of the generator as 
Optimum electric power generated 
7 Total input enthalpy flux 
Forming the ratio of (26) to the total input enthalpy flux, i.e. 
~ « 
m| ae jus], 
Y¥—-1po 
and transforming to dimensionless variables, we obtain 
O Rint (y/y —1) P(a) U(a) + 4U07%(a) zs 
peg ea a, (27) 
bk . Feoxt p 


where the last term represents the difference in entrance and exit stagnation 
enthalpy fluxes divided by the entrance flux. It is seen that the efficiency 7 is 
a complicated function (because of the complexity of the resistance ratio term) of 
the two dimensionless parameters 6 and M,, characterizing respectively the 
hydromagnetic and compressive properties of the system. 

Solutions of the system of equations (23) and (19) (evaluated at X = a) for k’ 
and U(a) were obtained on an I.B.M. 704 computer. Solutions were obtained for 
the particular Mach number M, = 0-3 and for a y equal to 3 (corresponding to the 
assumption that the plasma be treated as a monatomic perfect gas). d was allowed 
to vary continuously until a maximum dmax ~ 30 was reached beyond which no 
solutions could be generated. 
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Consider table 1 which lists the magnitudes of ali the dimensionless electrical 
and fluid mechanical quantities. It isseen that the solution exhibits an asymptotic 
behaviour near derit~ 4. In particular, if attention is focused on the column of 
exit Mach numbers (a), it is noted that for the second branch (d > 4), M(a) is 
always equal to one. We conclude, therefore, the existence of a finite range of 6, 
i.e. Ocrit < Od < dmax, where the required external loading at optimum power is that 
which forces a Mach number one flow at the exit. 


é - 2/k’ =« Rext/Rint I/Tine 1% U(a) P(a) M(a) 
0-1 20017 =0-99914 0-99166 1:0075 0-1466 1:0067 6:6096  0-30230 
0-2 2-0035 0-99823 0-98321 1-0153 0-29519 1-0136 6-5515 0-30468 
0-4 2-0076 0-99623 0-96585 1-0315 0-59858 1-0283 6-4321 0-30976 
0-6 2-0122 0:99394 0-94781 1-0487 0-$1076 1-044] 6-3080 0-31514 
0-8 2-0175 0-99131 0-92901  1-0671 1-2324 1-0612 6-1786 0-32101 
1-0 2-0236 0-98832 0-90942  1-0868 1-5642 1-0799 6-0434 0-32745 
2-0 2:0736 0-96450 0-79390 1-2149 3°4129 1-2093 5:2425 0-37296 
4-0 2:0077 0-99617 0-62218 1-6011 9-2909 2:7903 -6742 1-0000 


—_ 


6-0 1-4202 1-4083 1-2769 1-:1029 13-572 2-7237 1-6342 1-0000 

8-0 1-2407 1-6120 1-9160 0-84134 15-801 2-6882 1-6130 1-0000 
10-0 1-1562 1-7299 2-5444 0-67987 17-127 2-6671 1-6002 1-0000 
12-0 1-1083 1-8045 3-1645 0-57024 17-983 2-6533 1-5920 1-0000 
14-0 1-:0783 1-8547 3°7780 0-49093 18-565 2-6438 1-5863 1-0000 
16-0 1:0582 1-8899 4-3865 0-43085 18-976 2-6372 1-5823 1-0000 
18-0 1-0442 1-9154 4-9893 0-38389 19-274 2-6323 1-5794 1-0000 
20-0 1-0340 1-9342 5-5887 0-34608 19-496 2:6287 1-5772 1-0000 
25-0 1:0186 1-9636 77-0735 0-27760 19-845 2-6230 1-5738 1-0000 
30-0 1-0105 1-9791 8-5437 0-23165 20-358 2-6146 1-5687 1-0000 


TABLE 1. List of the dimensionless electrical and fluid dynamical variables. 
M,=0:3; y=% 


At least for the branch characterized by M(a) = 1, the condition which 
determines dmax may be formulated as follows. The Mach number at any point in 
the channel in terms of the dimensionless fluid dynamical variables P and U is 
given by M = ,(U/yP). The condition M(a) = 1 at the exit yields P(a) = U(a)/y. 
Substituting for P(a) in (23), we obtain 


U(a) _ (~+k'B)—[U (a) — ($k’) U%(a)) 
y 1—[y/(y—1)] k’U(a) 


é 


The simultaneous solution of (19), (23), and (28) yields dmax and the corre- 
sponding limiting values of k’ and U(a). 

Let us consider now the electrical quantities. Referring to the definitions of 
k and k’, both compressible and incompressible, it is seen that 2/k’ represents the 
ratio of the compressible terminal voltage to the terminal voltage obtained by 
assuming the gas to be incompressible. Similarly, J/Jinc, the ratio of the compres- 
sible to the incompressible current, is given by (2/k’) (Rint/Rext). It should be 
emphasized that incompressibility is defined and used here in the sense of § 4; 
i.e. the density of the fluid is everywhere constant in space and time and not as the 
limit flow approached as the Mach number is made to approach zero. 
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Plots of the resistance, voltage and current ratios as functions of the magneto- 
hydrodynamic interaction parameter 6 are shown in figure 2. Several general 
observations can be made. First, for weak interactions (i.e. d small), the fluid 
behaves as if it were incompressible; that is, the resistance, terminal voltage and 
current remain close to the corresponding incompressible values. Secondly, the 
asymptotic behaviour is clearly marked and is represented in the resistance and 
current curves as forming cusps at the critical 3. Thirdly, for large 6, the behaviour 
of the solution is very much different from the corresponding incompressible case. 
The terminal voltage ratio approaches the asymptotic value two while the external 
load ratio continues to increase sharply with 6. This behaviour is very significant, 
for it indicates (as is often done in fluid dynamics in order to obtain a first estimate 
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FicurE 2. Resistance, terminal voltage and current ratios against 
magnetohydrodynamic interaction parameter 6. M, = 0-3, y = 3. 


to some fluid mechanical interaction), that one cannot approximate the inter- 
action by simply treating the fluid as incompressible. For example, it is seen from 
the figure, that for 6 very large, the external load required at optimum power is 
close to an order of magnitude larger than the corresponding incompressible load. 

Figure 3 is a plot of the efficiency of the plasma generator as defined above. Asis 
expected, for small 5, when the effect of the magnetic field on the fluid flow is 
small, the efficiency is small. It increases as 6 increases approaching a value 
slightly over 20° in the neighbourhood of édmax. 

Finally, let us consider briefly the effect of the interaction on the fluid 
dynamical variables. The distribution of all of the fluid dynamical variables with 
distance along the channel was obtained for every 4 in each of the two branches 
with the view to determining whether any of the solutions belonging to the two 
branches corresponded to any peculiar phenomena, e.g. choking. In every case, 
the dynamical variables were single valued and varied continuously with channel 
distance so that behaviour peculiar to choking (Resler & Sears 1958) was absent. 
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Figure 4 shows the variation of the normalized fluid dynamical variables with 
fractional distance along the channel for the case 6 = 20. It is seen that the 
pressure and density decrease while the velocity and Mach number increases 
downstream of the entrance. It appears that, as far as the direction of change of 
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FIGURE 4. Variation of normalized fluid dynamical state variables with dimensionless 
distance along channel. M = 0-3, y = 3, 6 = 20. 


the fluid variables are concerned, the effect of the interaction is analogous both to 
channel flow with friction and no heat transfer or to channel flow with heat addi- 
tion and no frictional forces (see, for example, Shapiro 1953). The velocity increase 


along the channel leads to a very interesting conclusion concerning power delivery 
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to the external load. It was seen in § 3 that the contribution to the power delivered 
to the external load per unit channel length at station x is ko{uB—k/yo}. Since 
k, o and B are constants, the power contributed depends on the magnitude of the 
velocity at that station. We conclude, therefore, that most of the power delivered 
to the external load comes from the induction taking place near the exit region of 
the channel. 


6. Conclusion 

In this paper, solutions for the compressible generator were obtained for only 
one particular value of the entrance Mach number. It is planned to obtain 
a family of solutions over a wide range of Mach numbers both subsonic and 
supersonic. Interest here is not only toward the determination of the effect on the 
interaction of the various states of compressibility of the entering plasma, but to 
see whether the asymptotic behaviour occurs for other Mach numbers as well, and 
at what values of 6. The reason for, and the physical significance, if any, of the 
asymptotic behaviour itself should be investigated. Further, a thorough examina- 
tion should be made into the validity of the assumptions made in §1 directed 
towards estimating the practicality of the results. In particular, the effect of the 
possible development of space charge and the formation of space charge sheaths 
near the walls should be thoroughly investigated. 

In conclusion, it is hoped that, while this paper concerned itself mainly with the 
theoretical treatment of a very highly idealized but interesting problem in 
hydromagnetics, it will also serve in giving a first insight into the behaviour of 
a type of generator which may prove to have future practicality. 


The writer thanks Drs William McIlroy and Philip M. Mostov for their sugges- 
tions, clarifying discussions and helpful criticisms during the preparation of this 
paper. Thanks are also extended to Mr Sol Broder and Miss Eleanor Ditzler for 
their kind help in performing the numerical computations. 

This research was supported in part by the U.S. Air Force through the Air 
Force Office of Scientific Research. 


REFERENCES 

KuneNn, A. 1958 Electromagnetic pinch effect as a space propulsion system. Paper 
presented at the American Astronautical Society, Western Regional Meeting, Stanford 
University. 

NEURINGER, J. L. 1958 The formulation of a problem in optimum power generation using 
a plasma as the working fluid. Tech. Note Republic Aviation Corp. Sci. Res. Staff no. 17. 

Post, R. F. 1956 An application of the physics of high temperature plasmas. Rev. Mod. 
Phys. 28, 338. 

RESLER, E. L. & SEARS, W. R. 1958 The prospects for magneto-aerodynamics. J. Aero. 
Sci. 25, 235. 

SHaprro, A. H. 1953 The Dynamics and Thermodynamics of Compressible Fluid Flow, 
vol. 1. New York: Ronald. 





302 


Motion of a bore over a sloping beach 
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The results of numerical calculations are presented for the motion of a bore over 
a uniformly sloping beach. The shallow water equations are solved in finite 
difference form, and a technique is developed for fitting in the bore at each step. 
The results are compared with the approximate formula given by Whitham (1958) 
and close agreement is found. The approximate theory is considered further here; 
the main addition is a rigorous proof that, within the shallow water theory, the 
height of the bore always tends to zero at the shoreline. 


1. Introduction 


This paper gives an account of numerical solutions of the shallow water 
equations for the motion of a bore on a uniformly sloping beach. In a previous 
paper (Whitham 1958), an approximate formula was derived for the variation in 
the strength and height of the bore. The predictions were seen to contain specially 
interesting features, one of the main ones being that the height of the bore always 
tends to zero as it approaches the shoreline. Indeed, in the earlier account these 
predictions were viewed with a certain amount of suspicion, even though checks 
of the analogous theory in certain problems of gas dynamics show extremely high 
accuracy. Subsequently, it was realized that the predicted behaviour is correctt 
(within the shallow water theory at least) and the result that the bore height tends 
to zero at the shoreline can be proved rigorously. This proof, together with 
recapitulation and further discussion of the approximate method, is given in § 2. 

The numerical investigation confirms the previous results completely. In fact 
the agreement is quite remarkable; the differences are of the same order as the 
probable errors in the numerical work. In addition to this check, the numerical 
solution provides the full details of the flow behind the bore; these were not 
obtained previously. Furthermore, the solution is found for various starting 
conditions whereas it is assumed in the approximate theory that the bore is 
initially moving with constant speed and height in water of uniform depth. For 
these other starting conditions, it is found that the solutions eventually settle 
down and follow the approximate formula; the appropriate choice of the arbitrary 
constant factor in that formula varies considerably with the starting conditions 
for the same initial bore strength, but the functional dependence near the shoreline 
is the same. A similar insensitivity to the detailed initial conditions is well known 
for the analogous problem of the converging cylindrical shock in gas dynamics. 
The numerical results are discussed in detail in § 4. 


+ In this connexion we are grateful for discussions with G. F. Carrier and H. P. 
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For this problem, it was decided to develop a numerical method which uses fixed 
equal space intervals rather than the characteristics method which is not well 
suited to machine calculations. Stable ways of differencing the equations are well 
known, and the only question is how to fit in the bore at each step. This requires 
special care, but the method adopted here seems to be satisfactory in every 
respect. The details of the numerical scheme are given in §3. It is planned to 
extend this method for use in the analogous shock-wave problems in gas dynamics. 


2. Approximate formula and behaviour near the shoreline 

If h(x, t) = ho(x) +(x, t) denotes the depth of the water, where h,(z) is the 
undisturbed depth ahead of the bore, and u(x,t) is the particle velocity, the 
equations of the shallow water theory are 


m+ {(ho +7) Uh. = 9, (1) 
Um+UUu,+ gn, = Y. (2) 
The bore conditions are 

, h(hy +h) 
T — J? 0 ? 3 
2h (3) 

h—-h 

oar a (4) 


where U is the bore velocity. Along a positive characteristic it is readily found 
from (1) and (2) that 


TT 
du+2de—L = ¢ (5) 
U+C 
2 U+C 6 
on ° i 5 (6) 


where c = (gh). The approximate formula may be derived by applying the 
differential relation (5) to the flow quantities immediately behind the bore. Since 
these quantities are given in terms of U by (3) and (4) we have, then, a differential 
equation for U as a function of hj. The arbitrary constant in the integration is 
fixed from the initial value of U. 

This simple rule and its motivation are discussed in detail in the previous paper 
(Whitham 1958). It has not been justified in general, although it is easily shown 
to be the correct answer in a linearized perturbation theory for the effects of small 
changes in the depth on a bore initially moving with constant speed in a uniform 
region. In the linearized theory, the coefficient (wu +c)~! in (5) is replaced by its 
unperturbed value (uw, +c¢,)~? say. Then (5) can be integrated to 

Au + 2Ac— gAhy = 0, (7) 

U,+Cy 
where Au, Ac, Ah, denote the perturbations u — u,, etc.; the right-hand side of (7) 
is set equal to zero since all the positive characteristics C, come from the initial 
uniform region in which Au = Ac = Ah, = 0 (see figure 1). Since the combination 
in (7) vanishes on each characteristic, it follows that it vanishes everywhere. In 








304 H. B. Keller, D. A. Levine and G. B. Whitham 


particular, therefore, it applies to the flow quantities just behind the bore. 
Substituting from the bore conditions we have a relation of the form 


S(O) AU = Ahy (8) 


for the change in bore speed AU = U —U, in terms of the change in depth Ah. 
This result is exactly the linearized form of the above rule. The more general form 
of the rule corresponds to replacing (8) by the differential equation 
yee ii 
fie) = 1. (9) 
dhy 
From this point of view, the rule would be expected to hold if the depth varies 
sufficiently slowly even though the total change in depth is not small. But it turns 
out to be much better and applies equally well to quite extreme cases. This will be 
seen in the results presented in this paper. 





FicurE 1. The (x, t) diagram for the case u, < c¢,. 
1 1 


In substituting the bore conditions into (5) to get a differential equation 
corresponding to (9), it is convenient to work with MW = U/,/ (gh). In terms of MV, 


we have 


. 
= M .|(2M2-1), 
Ala ali 
u  —- 2M(M2-1) 
(gho) J (2.M?—1)’ 
NALA \ (10) 
, 
= ,/(2M?-1), 
Vigra) ! | 
| 
}) h —ho > | 
! = 2(M2~1), 
er (4° —1) 


0 0 


On substitution in (5), the following equation for Vas a function of hy is obtained: 


1 dh : (M+1)(M— 3)? (18+ M?- M —}3) - 
hodM- = (M—1)(M?—})(M4+3M3+ M?-3M-1) - 
Therangeof Wis 1 < M < «©, weak bores corresponding to small values of (.V — 1) 


and strong bores to large values of VW. In this range, the right-hand side of (11) is 





an 


W 





On. 


0) 


di: 





Motion of a bore over a sloping beach 305 


always negative so that M increases monotonically as h, decreases. However this 
is not true of the height 7. When M —1 is small, (11) is approximately 


Ldhp _ 4 1 
hdM 5 M—1" 


and it gives M-lechst, yohst. (12) 
Whereas, when J is large, we have 


| dh, 4 


a a Meahst. yochi. (13) 
_— a 


This is the result noted earlier that the bore height tends to zero when h, —> 0. For 
weak bores, 7 increases as h, decreases, but for strong bores 7 decreases with hyo. 
Thus in the beach problem, if 7 is initially less than the critical value 0-6262 hy (see 
the determination below), the height first increases and then ultimately decreases 
with a maximum at some intermediate position; at the maximum 7/h, = 0-6262. 
If 7/h) exceeds this value initially, the height decreases all the way to the shoreline. 
The bore velocity U has a similar change in behaviour, but at a different critical 
value given by 7/hg = 2-076, U/,/(gho) = 2-504. If U < 2-504, (gh) initially, 
U first decreases, but then ultimately increases to a finite value at the shoreline; 
at the minimum U/,/(gho) = 2°504. If U > 2-504,/(gho) initially, U increases all 
the way to the shore. 

The result that 7 > 0 as h, > 0 is surprising at first sight and certainly the 
validity of the rule might be suspected in this extreme limiting case. However, an 
independent proof can be given as follows: 

First, from the bore conditions (3) and (4), the bore height 7 must tend to zero 
with h, if w is to remain finite. We now show that u must remain finite. For 
a uniformly sloping beach we may take hy = «(x)—x), where « is the slope and 
x = X, is the shoreline. Then, the characteristic relations (5) and (6) give 


du+2dce+agdt=0 ie. w+2c+agt = constant, (14) 


on a positive characteristic. Thus, if «— 0 at the shoreline, u+2c+agt + a 
at all points of the characteristic which reaches the bore just at the shoreline. This 
is clearly impossible in the solution of our problem. (It should be noted that if 
uc, U +o so the bore reaches the shoreline in finite time.) Therefore, 
wremains finite as h, + 0; hence, from the bore conditions, 9 + 0. It is interesting 
to note that this method of fixing the behaviour near the shore, by the absence of 
a singularity on the characteristic, is completely analogous to the method for 
determining the exponent in Guderley’s similarity solution for the converging 
cylindrical shock (see the discussion in Whitham (1958)). 

Therelation (14) can be written ina further convenient form if weassume that the 
depth is uniform in z < 0and that the bore is initially moving with constant speed 
in that region. The (x, t) diagram is shown in figure | for the case when the initial 
constant values wu = u,, ¢ = c, behind the bore are such that wu, < c,. Then the 
flow is undisturbed to the left of the negative characteristic 2 = (wu, —c,)¢ shown 
as C_in the figure. If u, > c, the flow is undisturbed with u = u,,¢ = ¢,inz < 0 

20 Fluid Mech. 7 
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until a reflected bore appears. In the absence of such secondary bores (which 
would modify the result somewhat), we have from (14) 


u+2c+agt = u*(T) + 2c*(7T) + agr, 


where ¢ = 7 is the time when the characteristic crosses x = 0 and u*, c* are the 
corresponding values of u and c. But, w+ 2c is constant on the characteristic in 
x < Osince h, = constant, therefore u* + 2c* = u, + 2c, and we have 


u+2c+ag(t—T) = u,+ 2c. (15) 

ho = f(M) n l ¢ ue 

M A A V(gA) /(gA) V(gA) 
1-01 5-360 0-2155 2-385 2-361 0-092 
1-02 3-006 0:2429 1-839 1-802 0-137 
1-03 2-124 02586 1-590 1-543 0-173 
1-04 1-649 0-2692 1-440 1-385 0-202 
1-05 1-350 0-2767 1-339 1-275 0-228 
J 1-125 1-130 0-2825 1-261 1-188 0-252 
1-08 8-696 x 107 0-2894 1-163 1-077 0-290 
1-1 6-988 x 10-1 0-2935 1-096 0-996 0-324 
1-146 4-733 x 10-1 0:2964 1-005 0-877 0-387 
1:2 3°342 x 10-2 0-2941 0-951 0-793 0:445 
1:3 2-055 x 1071 0-2836 0-909 0-699 0-527 
1-428 1-285 x 1071 0: 2668 0-898 0-629 0-606 
1-5 1-028 x 10-3 0-2570 0-900 0-600 0-643 
2-0 3-308 x 10-2 0-1985 0-962 0-481 0°825 
/6 1-602 x 10-2 0-1602 1-028 0-420 0-935 
3-0 7-909 x 10-3 01265 1-100 0-367 0-035 
3°5 4-637 x 10-3 0-1043 1-155 0-330 1-106 
4-0 2-916 x 10-3 0-0875 1-203 0-301 1-164 
5-0 1-335 x 10-° 0-0641 1-279 0-256 1-253 
6-0 6-992 x 10-4 0-0489 1-337 0-223 1-318 
7-0 4-023 x 10-4 0-0386 1-383 0-198 1-368 
8-0 2-480 x 10-4 0-0312 1-420 0-177 1-408 
10-0 1-094 x 10-4 0:0217 1-476 0-148 1-468 
20-0 8-058 x 10-6 0-0064 1-605 0-080 1-603 
x0 0 0 1-763 0 1-763 


TaBLE 1. Calculations from approximate formula for bore motion: general A 


From this form it is immediately obvious that w remains bounded since all terms 
on the left are positive. However, this form is limited to the particular case of an 
initially uniform bore and the earlier argument is more general. With this result 
that w remains bounded, the bore conditions (3) and (4) show that h and 7 are 
proportional to he as h, -> 0. It should also be noted that the particle velocity u 
and the bore velocity U’ approach the same limiting value when hy > 0. 
Turning to the details of the solution given by (11), we see that the maximum 


of 7 occurs when 


dn 
dh) 


dM 7 


2(M2—1)+ A _= 
(VM 1)+ 4h, i. 


0. (16) 
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The value of is found to be 1-146, and from the bore conditions the corre- 
sponding 7, etc. are 
Maximum 9: M = 1-146, 9/hy = 0-6262, U/,/(gho) = it 
c , (17) 
~_ = 0-5627, = 1-275. 
v (gho) v(gho) 
(In the earlier paper (Whitham 1958), 1 was taken erroneously to be approxi- 
mately 1-2 and this led to a large error in 7/h,). The minimum bore speed occurs at 
a larger value of .V; we have: 


Minimum U: M = 1-428, y/hy = 2-076, U/,/(ghy) = 2-504, 


u c an (18) 
— = 1-690, —— = 1-754. 
v (gho) V (gho) 

H, = ho -_N= : V= . a v= ne ~ 

M hy (9) hy (0) V(gho(9)) V(gho(9)) J (gho(9)) 
J 1-125 1-0000 0-2500 1-186 1-118 0-237 
1-08 0-7696 0-2561 1-094 1-013 0-273 
EE 0-6184 0-2597 1-031 0-937 0-305 
1-146 0-4189 0-2623 0-946 0-825 0-364 
1-2 0-2957 0-2602 0-895 0-746 0-419 
1-3 0-1819 0-2510 0-855 0-658 0-496 
1-428 0-1137 0-2361 0-844 0-591 0-570 
1-5 0-0910 0:2275 0-846 0-564 0-605 
2:0 0-0293 0-1757 0-905 0-453 0-776 
4-0 0-0026 0-0774 1-131 0-283 1-095 
8-0 0-0002 0-0276 1-336 0-167 1-325 
00 0 0 1-659 0 1-659 

TABLE 2a. Calculations from approximate formula for bore motion: 
A = 0-8850, N(0) = 0-25 

Hy = Mo Newt se = —_ v= : 

M h,(9) ho(0) \ (gho(0)) \ (gho(9)) \ (gho(9)) 
/6 1-0000 10-000 8-12 3°317 7-39 
3-0 0-4936 7-897 8-69 2-897 8-18 
3°5 0-2894 6-511 9-13 2-608 8-74 
4-0 0-1820 5-459 9-50 2-375 9-19 
5:0 0-:0833 3:998 10-10 2-020 9-90 
6-0 0:0436 3°055 10-56 1-760 10-41 
7-0 0-0251 2-410 10-92 1-560 10-81 
8-0 0-0155 1-950 11-22 1-402 11-13 
10-0 0-0068 1-352 11-66 1-166 11-60 
20-0 0-:0005 0-401 12-68 0-634 12-66 
00 0 0 13-93 0 13-93 


TaBLE 2b. Calculations from approximate formula for bore motion: 
A = 62:41, N(0) = 10-0 


Equation (11) integrates to 
= (M? — 3) exp {a, tan (M +a,)/as} 
° (M=1)8 (M —a,)* (M? +a, M +a,)4 (M+a,)”" 
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where A is a constant of integration to be determined from the initial strength 


of the bore, and 


a, = 0-2808, a, = 0-6769, a, = 0-3179, a,=0-7471, a; = 1-354, 


a, = 05593, a, = 2393, a2 =1:180, 2 =1-173, y = 1-673. 


The flow quantities at the bore are determined from equations (10); their values 
are given in table 1. For a given strength of bore initially. MW is determined from 
the bore conditions and then A is determined from (19). The results for 9/hy = 0-25 
and 10 initially are given in tables 2a and 2b. They are discussed in detail in § 4. 


3. Numerical method 


It is assumed that the beach has slope a, starts at 2 = 0 and has the shoreline 
at x = 2); the depth is constant and equal to az, in x < 0. We introduce the 


dimensionless quantities: 





2 a 
X = XQ T= Xp * (ghy(0)). 
se u(x, t) r U 
) X ’ 1 = = ‘ 
me) = gh) V(gh(0)) 
- mm _ A(z, t) n(x, t) 
H(i ,7T) = h(0)’ N(L,T) = h(0)’ 
where h,(0) = aay. In terms of these variables 
HX) ={ a \ 
lA) = Da 0<X<1,J 


The shallow water theory equations (1), (2) can now be written as 
Hy, +(Hv)x = 9, 
Up + h(v?)y + Hy = Myx. 
and the bore conditions (3), (4) become 


H(H, +H) 
2H,’ 


7" (" Ih) - 


aur 


~ 


(20) 


The numerical procedure is, in outline, to compute the flow quantities behind 
the bore on a set of net points (X,, 7.) by means of finite difference approximations 
of (22) and (23). The bore quantities and its propagation between net points are 
computed from a special set of difference equations which include the bore 
conditions (24) and (25), and couple the flow variables at the bore to those behind 
it by means of the equations of motion. Of course this procedure assumes that 
the initial state contains a bore and in its present form is not capable of treating 


problems in which bores develop. 


In detail the spatial net is chosen to be uniform, X; = iA X, for convenience, and 
the time net, 7).,., = T,+A7Z,, must then satisfy the stability condition stated 
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below. If £,, is the position of the bore at time 7},, we let Xy,) be that net point for 
which " . 

X sx) S Sk S X 41 
At all net points X; < X,,) the flow quantities are computed from the difference 
forms of (22) and (23) (see figure 2a): 


AT 
H(P) = 3|H(R’) + H(Q’) — sry OR) A(R’) —0(@') HQ), (26) 
AT, pete Ra , P 
o(P) = dr(B’) + e(Q)) — ry A(R’) — (0)] + [H(R’) - HQ} 
0, X,; < 0, 
—14AT7,, X,; = 9, 
AT,, X;> 09. (27) 


x= &(T) 
lp Q P wo” 


Te +1 Te+1 








Ql P| R! Ql P! R! 


y 


























Xi-1 X; Xj+1 X;-1 X; Xse1 
FIGURE 2a. Mesh points for FiGuReE 2b. Mesh points for 
differencing the equations. fitting in the bore. 


These equations are obtained from (22), (23) by replacing the 7'-derivatives by 
forward differences and the X-derivatives by centred differences. If in the 
T'-differences the average values 3[H(R’) + H(Q’)], etc., are replaced by the actual 
values H(P’), ete., the difference equations become unconditionally unstable. 
However, using these average values it has been shown (Lax 1954) that the 
equations are stable provided the time mesh satisfies the so-called Courant 


condition | AY ) 


AT, < min ; ; 
. p \o(P’)+/{H(P')} 


/ 


It is clear from (21) that 0H,/¢X is discontinuous at X = 0 and as a consequence 
it follows that cv/0X and cH/cX must satisfy the jump conditions 


CH |X =+ —H =... v 
xl. =p lel. -z (28) 
E x-o— H-v E x-0o- H-v’ 


In these dimensionless variables, v? = H when the particle velocity is equal to the 
propagation speed. In fact this condition never occurs at x = 0; ifit holds initially, 
a reflected bore is produced when the main bore meets the change in depth, and 
this second bore establishes values of v and H, with v? + H, at x = 0. This whole 
question of perturbations on sonic conditions behind a shock has been elucidated 
by M. P. Friedman (1959). This particular feature of the problem is not under 
discussion here and we compute only cases with v?—H well away from zero. 
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Again the discontinuity in the bottom and corresponding discontinuities in vx 
and Hy, are not of primary interest here, since we are mainly concerned with the 
motion of the bore as it approaches the shoreline. Accordingly no special account 
of these discontinuities has been taken in (26) and (27). As a consequence we 
expect and indeed find in the results that the numerical solution suffers from 
relatively large truncation errors near X = 0. These effects could be eliminated 
either by smoothing the corner in the bottom profile or by modifying the difference 
equations to include the jump conditions (28). However, the errors introduced in 
this way at X = 0 are not large enough to warrant these detailed corrections for 
the present purposes. 

The net points which enter into the bore fitting procedure are shown in 
figure 2b. Allrelevant quantities are known at Q, Q’, P’ and #’ and the unknowns 
are o(P), H(P), v(R), H(R), V(R) and &,,, = €(R). Thus in order to have a deter- 
mined system we must add to the equations of motion and bore conditions two 
additional relations. We use the definition of the bore speed 





= V(T), (29) 


and the bore acceleration, obtained by differentiating (24) and using (22)—(25), 


dV(T) 1 


H H\? 
il od remade, THT ay : eae, ee? a pcos ‘ 
a 7 apt" (oH) x] (1+2 57) VAox (x) . (30) 


The equations of motion are approximated at P using backward T-differences and, 
for simplicity, X-differences centred at $(Q + R) rather than at P. These implicit 
equations are unconditionally stable and so impose no further restriction on the 
time mesh. The bore conditions (24) and (25) are to be satisfied at R and the 
subsidiary equations (29) and (30) are centred on the bore at }(R+R’). The 
system of equations obtained in this manner is non-linear and is solved by 
iterations. In the form and order in which these equations are used they are: 


E(R) = E(R’) + MAT, [V(R) + V(R)), (31) 
H(R) = H({Ho(R) [He(R) + 8V°(R)]} — Hy(B)}, (32) 
v(R) = 1 —ihn) V(R), (33) 
H(P) = H(P')~ ag) x UR) HR) —(Q) H@). (34) 
AP) = AP) gy cag — xy OR) — AQ] + LR) - HQ} 

0 £(R) < 0 

+ | wR x 4% E>0, X,< AX, 
AT, X, > AX, (35) 
AT, 
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Equations (32) and (33) are the bore conditions, equations (34) and (35) are the 
equations of motion and equations (31) and (36) are the difference forms of (29) 
and (30). Here W(R) and W(R’) are the lengthy but obvious difference forms of the 
right side of (30). In these evaluations the appropriate X-differences must be 
extrapolated to k and R’. By this procedure the fitting method yields the correct 
uniform flow over a uniform bottom. Any other procedure seems to require special 
considerations in order to obtain this steady solution. 

The iterations proceed from an initial estimate of V(R), say V(R) = V(R’), in 
(31). With the £(#) thus determined H,(R) = H,(&(R)) is computed from (21) and 
the quantities in (32)-(36) can then be evaluated in order. This procedure is 
repeated, using the latest value of V(R) to start each new iteration, until the 
successive iterates of all quantities are sufficiently close. In actual computations 
it is found that the number of iterations required to satisfy a prescribed con- 
vergence criterion decreases with the mesh widths and initial bore strength, and 
increases as the bore gets very close to the shoreline. 

When the iterations have converged the bore may have crossed one or more net 
points; i.e. €(R) > X,,, in figure 26. If this is the case, v and H at the intermediate 
net points are evaluated by linear interpolation between their values at P and R. 
To aid in smoothing the calculation the time mesh was restricted, in addition to 
the Courant condition, by AX 
In this manner at least two time steps are required for the bore to traverse a mesh 
width and the above interpolation procedure is used at most every other time 
step. For all of the calculations reported here it was observed that this latter 
condition determined A7’. 


4. Discussion of the results 


In the calculations a value is chosen for the initial ratio of the bore height to the 
undisturbed depth, 7/A,, and the other quantities at the bore then follow from the 
bore conditions. In the first calculations, the bore is started at x = 0 (or in some 
cases one or two mesh lengths to the left) with flow quantities behind the bore 
constant and equal to the values at the bore. The results are shown here for the 
cases with initial values 9/h, = 0-25 and 9/h, = 10. The first case is typical of 
those where 7/h, is less than the critical value 0-6262 so the height increases and 
the bore speed U decreases at first; also the initial flow behind the bore is subsonic, 
u, < ¢,. The second case is typical of strong bores with 7/hy above both critical 
values so that 7 decreases and U increases all the way to the beach; the initial 
flow behind the bore is supersonic, wu, > ¢,. 

The results for the dimensionless bore height 9/h (0) = NV and the dimensionless 
bore speed U/,/(gh,(0)) = V are shown in figures 3a and 3); V, indicates the value 
at the shoreline according to the approximate formula. In these calculations, the 
mesh width was AX = 0-01. Each calculation of the flow up to the arrival of the 
bore at the shoreline took about 45min on a Univac. It should be remembered 
that in the dimensionless variables the initial undisturbed depth h,(0) is 1, the 
shoreline is at X = 1 and the slope of the beach is 1. Any other depth h,(0), 
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distance x, or slope « (where ho(0) = a2), can be obtained from the scaling in (20), 
but the initial ratio of bore height to depth is fixed. In the first case, figure 3a, 
a small oscillation appears on the curve of V near £ = 0. This is the effect of the 
truncation error introduced by the discontinuities in the derivatives at X = 0 


r T T —————+ — —r — — 
| 
| 


0-30 | + 4 4 { hic ——{ 15 


0:20 


O10 





0 02 404 4206 O8 10 


FIGURE 3a. Variation of the bore height N and bore velocity V in the case N(0) = 0-25. 




















FIGURE 3b. Variation of the bore height N and bore velocity V in the case N(0) = 10-0. 


noted in §3. Corresponding disturbances in the other quantities could not be 
accurately detected in the numerical results. The observed oscillation in V has 
a magnitude of about 2 °% and is exceptionally large because of the small absolute 
value of NV in this case. As noted in §3 it did not seem worthwhile to modify the 
calculation specially to eliminate this particular inaccuracy. 

There is no point in plotting the corresponding values given by the approximate 
formula of § 2 since they lie almost on the curves. Instead the corresponding values 
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are given in tables 2a and 2b. Inthe case N(0) = 0-25, it is interesting to compare 
the predictions of the approximate formula for the maximum 7 and minimum U 
with the numerical calculations. They are shown in table 3. The numerical 
calculations can be taken with good accuracy to within a mesh width of the 
shoreline. In fact reasonable results are obtained at £ = 0-996. The limiting values 
of V at X = 1 may be obtained by extrapolation; this gives V = 1-568, for the 
case V(0) = 0-25 and V = 13-20, for the case N(0) = 10-0. Since V varies very 
rapidly in the last few mesh widths, we consider also an indirect method of 
extrapolation derived as follows: 

From the characteristic relation (14), we expect w+ 2c to vary little over the 
last few mesh lengths because the variation in ¢ is extremely small. Thus, since 
c > 0, a convenient prediction of uw at the shoreline is the value of w+ 2c (or its 


N _Y V - U 
E N Y Be & de J 
Maximum 7 Approximate formula 0-5811 0-2623 0-9457 0-6262 1-461 
Numerical calculation 0:5938 0-2587 0-9363 0-6367 1-469 
Minimum U Approximate formula 0-8863 0-2361 0-8444 2-076 2-504 
Numerical calculation 0-:8835 0:2349 0-8401 2-017 2-462 


TaBLE 3. Comparison of critical values 


extrapolation) over the last few mesh widths. As h, > 0, wu and U approach the 
same limiting value. The values obtained by this method for V at the shoreline are 
V = 1-681 for the first case, and V = 13-99 for the second case; the corresponding 
values from tables 2a and 20 are 1-659 and 13-93, respectively. In connexion with 
this extrapolation method it may be noted that the approximate characteristic 
rule also reduces approximately to u+ 2c = constant for strong bores. 

Turning now to details of the flow behind the bore, the height N and particle 
velocity v are shown for typical times in figures 4a and 4b. Again the effect of the 
truncation error introduced by the discontinuities at X = 0 is noticeable in 
figure 4a. In this first case, uw, < c, so that a reflected wave moves back into 
X < 0; in the second case, figure 4b, u, > ¢, so that the flow quantities in X < 0 
remain equal to their initial values. 

Figures 5a and 5b show the results of different starting conditions. The initial 
bore height is kept at (0) = 0-25, as for the case shown in figure 3a, but behind 
the bore two different initial distributions of V and v are chosen in place of the 
constant values. First NV and v are taken to fall linearly in —0-1 < X < Oto {4 of 
their values at the bore, and then remain constant in X < —0-1; secondly, NV and 
v fall in the same distance to } their values at the bore. The results for the bore 
height N and bore speed V are shown as curves 2 and 3 in figures 5a and 5b, 
respectively, and compared with the original curves | of figure 3a. At the points 
where V is minimum the values of V/,/H, = U/,/(gho) are 2-462, 2-539, and 2-674, 
for curves 1, 2 and 3, respectively. These are reasonably close to each other and to 
the value 2-504 given by the approximate formula in (18). This indicates that the 
formula (19) still applies near the shore even though the starting conditions are 
not the ones assumed in that approximate theory. Of course the values of A will 
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be different for the three different cases shown. Curves 1 have A determined from 
the initial bore height according to §2. For curves 2 and 3, the approximate 
formula does not apply in the early stages and some other choice of A must be 
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FIGURE 4a. The flow behind the bore: height N and particle velocity v 
in the case N(0) = 0-25. 
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Ficure 4b. The flow behind the bore: height N and particle velocity v 
v in the case N(0) = 10-0. 


made to see if the formula can describe the curves near the shore. It is convenient 
to choose A in each case so that the minimum of V agrees. The corresponding 
curves are drawn dashed in figures 5a and 5b using the values in table 1 with the 
appropriate choice of A, which is 0-8004 for curves 2 and 0-5325 for curves 3. 
These show to what extent the motion of the bore depends on the detailed initial 
distribution of NV and V, before tending to one of the curves (19) as the shoreline is 
approached. It is interesting to see (figure 5a) that in each case the curve goes 
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approximately to the value in the uniform region behind the linear profile. This 
might have been expected of course. It should be remarked that these initial 
constant values of NV and v in X < —0-1 do not satisfy the bore conditions 
exactly; for VN = 0-225, the bore conditions give v = 0-2627, whereas it was taken 
to be 0-2135; for N = 0-125, the bore conditions give v = 0-1367, whereas it was 
taken to be 0-1186. 




















FicureE 5a. Effect of different starting conditions: bore height NV. 
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Ficure 5b. Effect of different starting conditions: bore velocity V. 


Finally figures 6 and 7 show on a greatly magnified scale, the oscillations due to 
the truncation errors at X = 0 for the case N(0) = 0-25. These calculations were 
run until the bore reached £ = 0-159 for different mesh sizes and the convergence 
as AX - Ois clearly indicated. It was observed (figure 7) that while the mesh size 
changed by a factor of 8 the value of N at the bore changed by only 4 %. The value 
AX = 0-01 was used for the full calculations described above. 
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All these calculations were performed on the A.E.C. Univae at New York 
University. 


This research was supported by the United States Air Force under Contract 
No. AF 49(638)-446, monitored by the AF Office of Scientific Research of the Air 
Xesearch and Development Command. 
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FicuRE 7. Convergence of oscillation in height V(X, ¢). 
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REVIEWS 


Applied Hydrodynamics. By H. R. Vavtentine. London: Butterworth’s 
Scientific Publications, 1959. 272 pp. 50s. 

In the preface the author states that the book is intended ‘as an introduction to 
Hydrodynamics for students of Applied Mathematics and as a course in Fluid 
Dynamics for senior and post graduate students in Civil, Mechanical and Aero- 
nautical Engineering’. He also states that ‘the book has a bias towards practical 
application’ and that ‘particular attention has been paid to the provision of 
detailed physical explanations’. 

Broadly speaking, it can be said that the author has written a book which is 
suitable for senior engineering and applied mathematics students as regards 
level and treatment, but there are serious omissions in the subject matter. Each 
aspect is taken slowly and illustrated by worked examples. From a mathe- 
matical stand-point the treatment is reasonably rigorous and the physical back- 
ground is kept well in mind. References are given to further reading and there 
is a list of problems (with answers) at the end of each chapter. The book is 
limited to incompressible flow and, apart from a short chapter on the effects 
of viscosity, is mainly concerned with the irrotational flow of an inviscid fluid 
in two dimensions. It is attractively produced, being well illustrated with 
diagrams, and the use of headings in heavy type makes it easy to find any 
particular subject. 

Chapter 1, headed ‘The flow of an ideal fluid’, introduces a variety of topics, 
including fluid properties, continuity, stream lines, streak lines and path lines 
(or filament lines), rotational and irrotational flow, the stream function and 
velocity potential, Laplace’s equation, flow nets, Euler’s equation for inviscid 
flow, Bernoulli’s equation and its applications. The treatment of these topics is 
in keeping with the author’s intentions and is generally satisfactory, but Kelvin’s 
theorem (relating to the circulation in a circuit moving with the fluid) should 
have been proved after the derivation of Euler’s equation for inviscid flow. There 
is only a passing reference to this important theorem in a later chapter. 

Chapter 2 is entitled ‘The flow of a real fluid’ and deals with the effects of 
viscosity, with the distinction between laminar and turbulent flow, and with 
boundary layers. Results for velocity distributions are quoted and not derived, 
which is unconvincing for the student. I consider this the least satisfactory 
chapter in the book, largely because too little space is devoted to it, with the 
result that the discussion of the various themes is inadequate. For instance, in 
dealing with stresses in a viscous fluid there is no reference to the deformation of 
a fluid element and the rates of strain. A fuller discussion of the laminar and 
turbulent boundary-layer flows, from both theoretical and experimental aspects, 
would have helped this chapter very much, and the derivation of the two- 
dimensional equation of motion might have been given in an appendix. These 
modifications would require a good deal of space and it might be argued that 
such a treatment is outside the scope of the title of the book, but the book is 
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intended for engineers who have to deal with real fluids; in any case, where a 
problem can be treated by inviscid flow methods it is important to understand 
why this is so! 

Chapter 3, on ‘Graphical flow nets, numerical analysis and experimental 
analogies ’, is a very welcome chapter which is quite well done. There is a variety 
of examples of the appropriate construction of flow nets by graphical methods and 
a short account of the solution of Laplace’s equation by numerical methods. There 
is a useful account of the membrane analogy and of Hele-Shaw’s viscous flow 
analogy, but more space should have been given to the electrical analogy. This is 
the most useful analogy for engineers, since it is capable of yielding fairly precise 
quantitative information and it can be applied to a wide variety of problems. 

Chapter 4 deals with the ‘Standard patterns of flow’. There are both graphical 
and analytical treatments of the elementary flows such as source, sink, doublet, 
point vortex and the combinations with a uniform stream. The concept of 
circulation, which is introduced here for the first time in the book, would have 
been much better placed in the first chapter of the book. 

Chapters 5 and 6 are headed ‘Conformal transformation’ and together occupy 
about one-third of the book. Chapter 5 introduces complex numbers and 
functions of a complex variable leading to the solution of Laplace’s equation in 
two dimensions. There is a reasonable account of analytic functions and of 
singular points, but, surprisingly, nothing is said about multi-valued functions, 
the need to ensure a point-to-point correspondence between corresponding 
regions of the various planes, and the use of ‘cuts’ or Riemann surfaces. This I 
regard as a grave defect, since in my experience this is where engineering 
students often have difficulty. The potential functions for the elementary flows 
are established in this chapter, and, finally, the Joukowski transformation is used 
to derive the flow about a normal flat plate and the Joukowski family of aero- 
foils. The chapter ends with an account of the lift on a two-dimensional aerofoil. 
I think it is a pity that the pressure distribution round an inclined flat plate at a 
small angle of attack was not dealt with in detail, because this is fundamental to 
the understanding of flow past an aerofoil. Chapter 6 is devoted to flows in- 
volving straight boundaries. The Schwarz—Christoffel theorem is stated and its 
meaning is explained. It is then applied to the mapping of semi-infinite and 
infinite strips and a number of flows involving straight boundaries are dealt 
with in detail. These include sink flows in infinite and semi-infinite channels, the 
flow past a normal flat plate and the flow into a rectangular channel. Finally, 
an account is given of free streamline theory and there are examples of its 
application, including the flow through a slot in a plane wall, Borda’s mouth- 
piece and flow past a normal flat plate. 

The last chapter is headed ‘Three-dimensional irrotational flow’. It is really 
devoted to axisymmetric flows, and it begins by developing the equation of 


continuity and the equations governing the stream function and velocity 
potential for cylindrical and spherical polar co-ordinates. The elementary flows 
are then developed and followed by their combination with the flow due to a 
uniform stream. Lastly, there is a short account of the numerical analysis of 
axisymmetric flows. 
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At the end of the book there are three appendices, one on the theorems of 
Green, Stokes, Cauchy and Blasius (to which Kelvin’s theorem might have been 
added), one on exponential, trigonometric and hyperbolic functions (useful 
although hardly necessary), and one summarizing the main relations and 
equations governing the flow of an inviscid incompressible fluid. 

In addition to the omissions I have mentioned under various chapter headings, 
there are two others which I regard as very important in a book at this level 
intended for engineers. First, I should like to have seen some account of the 
approximate methods which are used in aerofoil theory; for example, the 
development of fairings using continuous distributions of sources and sinks, and 
the development of camber line theory using the concept of a bound vortex 
sheet. Secondly, there should have been an account of vortices in three dimen- 
sions, together with their application to the theory of wings of finite span. 

My conclusions are that as a text-book on applied fluid dynamics it is too 
narrow in the selection of subjects, particularly for mechanical and aeronautical 
engineering students. A number of practically important aspects of inviscid 
incompressible hydrodynamics have not been dealt with and the effects of 
viscosity merit a more detailed and comprehensive treatment. On the other 
hand, there are useful aspects of the subject which have been dealt with and 
which are not usually found in text-books. The level and treatment of the 
subjects selected are entirely suited to the students which the author had in 
mind and there is much in the book which will be of great value to them, even 
though it cannot be regarded as a comprehensive text-book for a course in 
incompressible fluid dynamics. J. H. Panerox 


The Potential Theory of Unsteady Supersonic Flow. By J. W. Mizzs. 
Cambridge University Press, 1959. 220 pp. 45s. 
This book is the eleventh in the series of Cambridge Monographs on Mechanics 
and Applied Mathematics (edited by G. K. Batchelor and H. Bondi) to be issued, 
and is in a sense complementary to the monograph by Ward on Linearised 
Theory of High-Speed Flow which was published some time ago. Whereas, 
however, many books have been written on high-speed steady flow, no compre- 
hensive treatise exists on the theory of unsteady flow. For this reason this 
monograph will be welcomed by all who are concerned with flutter problems 
which require a knowledge of unsteady aerodynamic forces at supersonic speeds. 
It supplements the general information on unsteady flow problems included by 
I. E. Garrick in Section F, entitled ‘Nonsteady wing characteristics’, of Vol. vir 
of High Speed Aerodynamics and Jet Propulsion, by A. T. van de Vooren in an 
article on ‘Unsteady aerofoil theory’ published in Vol. v of Advances in Applied 
Mechanics, and by G. Temple in a chapter written for Modern Developments in 
Fluid Dynamics: High Speed Flow. Related information is also to be found in 
comprehensive treatises on aeroelasticity by R. L. Bisplinghoff, H. Ashley and 
R. L. Halfman and by Y. C. Fung. In addition, a recent series of papers on 
unsteady transonic flow theory by M. Landahl contains much valuable informa- 
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tion. There still remains, however, a need for a comprehensive account of 
unsteady subsonic theory. 

The present monograph is mainly concerned with supersonic flow and only 
draws attention to corresponding subsonic solutions where comparison is of 
value in giving further information on the problem. It starts with a discussion 
of the basic assumptions made in the development of the linearized equations 
used in the theory, and the validity of various approximations are considered 
by extending the method of analysis used by Lin, Reissner and Tsien in their 
classic paper on two-dimensional non-steady motion of a slender body in a com- 
pressible fluid (J. Math. Phys., 1948). Then follows a discussion of various trans- 
formations which have been used in subsonic theory and it is shown how similar 
methods can be used in the theory of unsteady supersonic flow. Supersonic 
problems are much simpler than those of subsonic flow because the effect of the 
wake can be neglected when the component of flow normal to the trailing edge 
is supersonic. In many cases exact solutions can be obtained, and four chapters 
are devoted to a discussion of such solutions. For many of the problems alter- 
native solutions exist but most of those described are due to the author who has 
made very many important contributions in this field of research. 

The chapter on slender wing methods and their application to unsteady flow 
problems is fairly comprehensive and is a most valuable part of the monograph. 
In the following chapters solutions are obtained for oscillating delta wings, low 
aspect ratio rectangular wings and slender non-planar bodies. The monograph 
concludes with a brief discussion of non-linear problems, and reverse flow theorems 
are discussed in the Appendix. A very extensive and useful list of references is 
included. 

The monograph is an authoritative account of the subject and one that should 
be studied by all interested in unsteady supersonic flow theory. 

W. P. JoNEs 


Proceedings of the Sixth Midwestern Conference on Fluid Mechanics. 
The University of Texas, 1959. 465 pp. $12.50. 


The above conference and the Fourth Midwestern Conference on Solid Mechanics 
were held simultaneously at the University of Texas in September 1959, and 
copies of both sets of proceedings are now available. The original typescripts 
and drawings submitted by authors have been reproduced photographically, 
which certainly makes for speed of publication although not for attractiveness 
of appearance. The volume on fluid mechanics contains 29 contributed papers 
and the text of an invited lecture by S. F. Schaaf on rarefied gas dynamics. 
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